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ABSTRACT 

^ . This is the first paper of a series that describes the methods and basic results of the 

\Q ' GALICS model (for Galaxies In Cosmological Simulations). GALICS is a hybrid model for 

OO ! hierarchical galaxy formation studies, combining the outputs of large cosmological N- 

t-H ■ body simulations with simple, semi-analytic recipes to describe the fate of the baryons 

ON ' within dark matter halos. The simulations produce a detailed merging tree for the 

dark matter halos including complete knowledge of the statistical properties arising 
from the gravitational forces. We intend to predict the overall statistical properties 
of galaxies, with special emphasis on the panchromatic spectral energy distribution 
. emitted by galaxies in the UV/optical and IR/submm wavelength ranges. 

O i 1 In this paper, we outline the physically motivated assumptions and key free pa- 

rameters that go into the model, comparing and contrasting with other parallel efforts. 
We specifically illustrate the success of the model in comparison to several datasets, 
showing how it is able to predict the galaxy disc sizes, colours, luminosity functions 
from the ultraviolet to far infrared, the Tully-Fisher and Faber-Jackson relations, and 
the fundamental plane in the local universe. We also identify certain areas where the 
model fails, or where the assumptions needed to succeed are at odds with observations, 
and pay special attention to understanding the effects of the finite resolution of the 
'jj] ■ simulations on the predictions made. Other papers in this series will take advantage of 

different data sets available in the literature to extend the study of the limitations and 
predictive power of GALICS, with particular emphasis put on high-redshift galaxies. 
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1 INTRODUCTION 

In the last five years, the discovery of the Cosmic Infrared 
Background ^Puget et jil^2£96j^^|C_u jderdoni et al. 19971 
IFixsen et al. 19981 Mauser et al. 19981 and the faint 
galaxy counts with ISO at 15 fim llElbaz et al. 19991 
and 175 /xm (|Puget et al. 1999) , SCUBA at 850 /im 
Smail, Ivison, fc Blain 1997| and MAMBO at 1.3 mm 
Carilli et al. 20001 have shown that about two-thirds of 
the luminosity budget of galaxies is emitted by dust in 
the IR/submm range. Whilst the nature of the source 
that heats up the dust is still uncertain, it now seems 
increasingly plausible that the contribution of AGNs to 
heating is not dominant, and that most of the energy 
is powered by starbursts due to gas inflows triggered by 
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close encounters and merging. The IR/submm wavelength 
range is actually tracking the star formation rate history 
of the Universe more accurately than the UV range, with 
strong sensitivity to the merging phenomenon that is the 
signpost of hierarchical galaxy formation. The luminous 
and ultra-luminous infrared galaxies that contribute to 
the CIB are thought to be the progenitors of the bulges 
and elliptical galaxies in the local Universe. The goal of 
GALICS is to get a consistent panchromatic description of 
this process and of the luminosity budget of galaxies as it 
appears, for instance, through the faint galaxy counts at 
various wavelengths in the optical, IR and submm. 

Various pieces of work have converged to build up a con- 
sistent description of galaxy formation within the paradigm 
of hierarchical clustering. Initial density perturbations 
are gravitationally amplified and collapse to form almost 
relaxed, virialized structures called dark matter halos. In all 
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variants of the Cold Dark Matter scenario (IPeebies 19821 
IBlumenthal et al. 19 84). smaller halos form first, and bigger 
halos form continuously from the collapse of smaller halos. 
Gas radiates and cools down in the potential wells of the ha- 
los (TWhite fc R ecs 19781. Halos have little angular momen- 
tum, and dissipative collapse stops when the cold gas settles 
in rotationally-supported discs (IFall fe E fstathiou 19801 
Dalcanton, Sperge Tj" fc Summers 1997 
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scales can be reasonably described with Schmidt laws 
or more sophisticated recipes (Kennicutt 1989). Models 
of spectrophotometric evolution of the stellar popula- 
tions produce luminosities, spectra, and colours from the 
star formation rate histories and initial mass function 
ijBruzual 19811 |Guiderdoni fc Rocca-Volmerange 1987| 
IBruzual fc Chariot 19931 . and are now extended 
to the modelling of dust thermal emission 

l|Mazzei, Xu, fc de Zotti 1992| ISilva et al. "19981 

|Devriendt, Guiderdoni, fc Sadat 19991 hereafter DGS). 
When they die, stars eject gas, heavy elements and 
energy into their environment. The energy feedback 
heats up the remaining gas and can produce galac- 
tic winds in the shallower potential wells that deplete 
the galaxies and quench subsequent star formation 
lIDekel fc Silk 19 86 1 . Finally, spheroids form from ma- 
jor mergers lIToomre fc Toomre 19721 IKent 1981) . If the 
spheroid can still be a centre for gas cooling, a new disc 
forms around this bulge. As a result, the morphological 
and spectral types of galaxies are not fixed once for all, but 
rather evolve as star formation, gas accretion and merging 
occur. 

These ingredients can be put together with some 
success within a fully semi-analytic model (hereafter SAM) 
that starts from the power spectrum of linear fluctuations, 
and follows the various processes right up to spectral energy 
distributions of stellar populations. For instance, under 
the assumption that each newly-collapsed peak produces 
a halo where a new galaxy forms, therefore neglecting 
the classical 'cloud-in-cloud' problem, it is possible to 
reproduce at least qualitatively the main statistical prop- 
erties of galaxies ijLacey fc Silk ^^ij ^ |Lacey et al. 1993| 



IGuiderdoni et al. 19981 fHeroendT fc Guiderdoni 200011 
The Extended Press-Schechter prescription (EPS, 
IPress fc Schechter 19741 IBond et al. 19911 IBower 19911 
|Lacey fc Cole 1993) is a more efficient tool to describe 
gravitational collapse and estimate the merging' history 
trees of dark matter halos. The fate of the stars, gas and 
heavy elements can also be followed within the hierar- 
chy of merging halos with an implementation of these 
ingredients in the EPS formalism ([White fc Frenk 199111 . 
However it is only through Monte-Carlo realizations of 
halo merging history trees (JRauffmann & Whit e~1993l 
|Lacey fc Cole 1994| ISomerville fc Kolatt 19991 that galaxy 
merging can be followed and hierarchical galaxy formation 
can be addressed. Implementing a simple recipe for dy- 
namical friction of satellite galaxies in the potential wells 
of ha los enabled [Kauffmann, White, fc Guiderdoni ( 1993| > 
and |Kauffmann, Guiderdoni, fc White (1994{ to follow 
the galaxy merging history. Though the 'block model' for 
hierarchical structure formation has been used for some 
time HCole et al. 19941 . most studies now involve this type 
of random realization based on the EPS for the merging 



history (|Somerville fc Primack ( 19991, hereafter SP99; 
|Cole et al. (2000) , hereafter CLBF, and papers of these 
series). In addition to the cosmological parameters (Ho, fio, 
Qa, f2s, the shape of the power spectrum and its normal- 
ization erg), the semi-analytic method introduces a limited 
set of free parameters because some of the processes have to 
be addressed phenomenologically: in the most general case 
they can be reduced to a star formation efficiency, a stellar 
feedback efficiency for the ISM/IGM, and a parameter 
that describes our ignorance on the complicated merging 
processes. These parameters are determined by requiring 
the results to fit certain datasets, commonly including 
the A"-band luminosity density in the local Universe, the 
number of dwarf galaxies (which are the most sensitive 
objects to stellar feedback) and the number of elliptical 
galaxies (which only form from major mergers in the 
simplest hierarchical scenario). Once the free parameters 
are fixed, many predictions can be produced and compared 
to data. 

A large number of papers have been devoted to var- 
ious aspects of galaxy formation and evolution, ranging 
from the Butcher-Oemler effect iKauffman n 19951 . the for- 
mation of discs and bulges (IKauffmann 19 96a I. damped 
Lyman-Q systems (Kauffmann 1996b), Lyman- break galax- 
ies jBaugli et al. 1998"l|Somerville, Primack, fc Faber 2001) , 
and the parallel evolutions of quasars and galaxies 
UKauffmann fc Haehnelt 200011 . 

However, this approach still suffers from a number of 
shortcomings. First, even if the EPS agrees with Af-body 
simulations ( Efstathio u et al. 19881 
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Kauffma nn fc White 19931 
ISomerville et al. 2000) . it is clearly a limiting simplifi- 
cation of the complex dynamical processes that actually 
occur. The non-linear dynamics is computed with the 
top-hat model that assumes sphericity and homogeneity, 
and it overestimates the number of halos on galactic and 
group scales ijOross et al. 1998|l . The only pieces of spatial 
and dynamical information that are stored in the halo 
merging history trees are the virial radii and circular ve- 
locities. There is no information on the spatial distribution 
and peculiar velocities of halos. As a consequence, the 
outputs of SAM cannot be used for synthesizing realistic 
mock catalogues and images that take into account spatial 
correlations, whereas there is an increasing need for these 
catalogues and images to analyse current and forthcoming 
observations and to test data processing techniques. 

It is tantalizing to bypass some of these limits by us- 
ing merging history trees produced from large cosmological 
A^-body simulations. The basic idea is to get a description 
of the dark matter halo merging history trees, which can 
be computed accurately from simulations, and to keep the 
usual semi-analytic approach to model the more uncertain 
physics of baryons. This model rests on the assumption that 
baryons do not alter significantly the dynamics of the dark 
matter, except on the smallest scales. Hence the name 'hy- 
brid model'. The result is a more realistic merging history 
for halos, which necessarily reflects on galaxy formation and 
evolution. The drawback is a loss of flexibility since the val- 
ues of the cosmological parameters Ho, fioj and £7 a, as well 
as the choice of the power spectrum P(k) and normalization 
erg, are built in the Af-body simulation. However, the value 
Qb, the physics of baryons, and the associated free param- 
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eters can be changed ad libitum, very much as in classical 
semi- analytic models. 

The first attempts at a hybrid model were proposed by 
IWhite et al. 19871 and Roukcma ct al. 1997 Less than ten 
snapshots of A/-body simulations were considered at that 
time, and this crude time resolution had a heavy impact 
on the results. Moreover, it is also obvious that any numer- 
ical simulation has a finite mass resolution, and that the 
unknown fate of baryons in systems below this threshold is 
going to propagate over the threshold in the picture of hier- 
archical clustering. A partial solution of this problem is to 
use the spatial information of a numerical simulation, but 
to build halo merging history trees from Monte-Carlo real- 
izations of the EPS i Kauffmann, Nusser, fc Steinmetz 1997| 
IBenson et al. 20001 Govern^o^t^iL^00Ty ~T Jnfortunatelv. 
it becomes impossible to follow the evolution of galaxies 
backwards in the structures, and the merging trees are 
not consistent with the merging histories of the individ- 
ual halos in the simulation. Only fully hybrid models keep 
this record and have fully-consistent merging history trees 
( Kauffma nn et al. 19991 hereafter KCDW, and following pa- 
pers). Moreover, only most recent hybrid models, taking ad- 
vantage of very high resolution TV-body simulations to dy- 
namically follow substructure within dark matter halos, are 
capable of accurately tracking properties of individual galax- 
ies within clusters jSpringel et al. 2001| > . 

Here we propose the GALICS model of hierarchical 
galaxy formation which is intended to provide a fully 
panchromatic description of the galaxy merging history, sim- 
ilar in spirit to that implemented by Grana to et al. 20001 in 
a pure semi-analytic model. For that purpose, we will follow 
chemical and spectrophotometric evolution in a consistent 
way, estimate dust extinction and radiation transfer, and 
compute spectral energy distributions of dust thermal emis- 
sion, following the lines of the STARDUST model (DGS). 

Our main goals are to: 

• present an original hybrid model that is entirely inde- 
pendent of previous attempts, and study its successes and 
failures compared to other models. 

• present an overall view of galaxy evolution by producing 
a host of predictions from a 'standard' reference model to 
compare with observed galaxy properties locally and at high 
redshift. 

• produce a panchromatic picture, which is closely linked 
to the hierarchical formation of structure, as galaxy mergers 
are bright in the infrared. 

• implement the effects of observational selection criteria 
that follow as closely as possible the actual observational 
processes. 

• study the effect of mass and time resolution constraints 
on these predictions. 

This paper (the first in a series) proposes an overall pre- 
sentation of the model and the basic predictions for those 
statistical properties of galaxies in the local Universe that 
can be more easily compared with other works. In section|5| 
we describe the procedure that has been developed to find 
the halos and build the halo merging history trees from the 
series of output snapshots. Section [3] presents the 'recipes' 
for cooling of hot gas in the halos. In section 01 we intro- 
duce a simple, standard implementation of the dissipative 
physics of baryons, and the construction of the galaxy merg- 



ing history trees. Section |K| describes the modelling of merg- 
ing events, and how these events are assumed to drive galaxy 
morphologies. In section|H]we present our methods for deter- 
mining galaxy luminosities, which are largely based on the 
STARDUST model of DGS. Section briefly summarizes the 
few free parameters that go into the models. In section|H] wc 
show results for the z = properties of galaxies in our ref- 
erence model, including sizes, colours, luminosity functions, 
and structural relations (Tully-Fisher, Fab er- Jackson, Fun- 
damental plane). In section [5] we investigate in detail the 
effect that the finite spatial, mass and time resolution of 
our simulations has on the predictions of galaxy properties. 
Section 1101 presents a discussion of these first results. Ap- 
pendix^describes the parallelized treecode that we use for 
large cosmological A/-body simulations. 

Four other papers will complete this description of the 
model. In paper II (Devriendt et al. in preparation), we will 
explore the sensitivity of our reference model to changes 
in some of the modelling assumptions, and study the in- 
fluence of variations in the choice of cosmology, recipes for 
the baryonic processes, and astrophysical free parameters. 
This paper will also discuss the evolution of galaxy prop- 
erties. In paper III (Blaizot et al. in preparation) we will 
focus more specifically on predicted properties for Lyman- 
break galaxies at redshift 3. In paper IV (Devriendt et al. in 
preparation), we will give faint galaxy counts, source counts, 
and angular correlation functions in the UV, optical, IR and 
submm ranges. Finally, in paper V (Blaizot et al. in prepara- 
tion), we will focus on the redshift distributions of different 
classes of galaxy, the redshift evolution of the spatial correla- 
tion function, £(r), and statistical bias, b, and the influence 
of environment on galaxy properties. Forthcoming papers of 
the series will address several other issues with an improved 
modelling of the baryonic processes. 



2 DARK MATTER HALOS 

We have used for the simulation a parallel tree-code writ- 
ten and optimized for the Cray T3E. This code, based on 
hierarchical methods, is described in detail in |Ninin (1999| l. 
We assume that, initially, the baryon density field traces 
that of the dark matter, and thus that galaxy formation 
occurs at local maxima of the underlying density distribu- 
tion. For a given cosmology it can be shown that there is 
a certain turn-around density contrast, above which mat- 
ter has collapsed into gravitationally bound systems which 
have separated from the expansion of the Universe. The pre- 
cise density contrast at which to define a 'virialized' halo is a 
matter of some debate i White 2001L and in general depends 
on cosmology. A physically sensible definition is that within 
the virial radius of an object the dark matter is virialized, 
and external to it material is infalling. This is found to occur 
at around 200 times the critical density, regardless of cos- 
mology. It is this latter definition we will use in this work. 
We refer to regions attaining this overdensity as dark matter 
halos, and it is assumed that all galaxy formation processes 
take place within these halos, and furthermore that there is 
no communication between these halos. 
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150.0 


Particle mass[lO 11 M0] 


0.08272 


Omega matter f2o 


0.333 


Omega lambda Hyv 


0.667 


Omega curvature O c 


0.0 


Hubble parameter h 


0.667 


Variance ag 


0.88 


r = Q ^ 


0.22 


Initial redshift 


35.59 


Number of timesteps 


19269 


Number of outputs 


100 


Spatial resolution (kpc) 


29.29 


Work time (10 3 h) 


56. 



Table 1. Simulation parameters used in our standard (ACDM, 
Apart = 256 3 ) simulation. 



2.1 Simulation details 

We outline the details of the treecode approach in ap- 
pendix^ We have initially simulated a spatially flat, low- 
density model (ACDM) with a cosmological constant, in a 
cube of comoving size Lbox = 100/i -1 Mpc. The amplitude 
of the power spectrum was set by demanding an approx- 
imate agreement with the present day abundance of rich 
clusters ( |Eke, Cole, fe Frenk 1996|l. Initial conditi ons were 
obtained with the GRAFIC code flBertschinger 1995| > . A loga- 
rithmic spacing in the expansion parameter, a, was used for 
the output file times: we finally had around 100 output files. 

The simulation contains A'part = 256 particles. We 
ran additional simulations at lower resolution (A'part = 
128 3 ,64 3 ), which we will use to test for convergence in the 
properties of our dark matter halos and galaxies (see sec- 
tion EJ. The full details of our standard simulation are sum- 
marised in table Q 



2.2 Detection of the halos in the simulation 

Many algorithms exist for detecting halos of dark matter in 
Af-body simulations. Examples include denmax, the spheri- 
cal overdensity algorithm, and the 'hop' algorithm. None of 
them, in spite of their complexity, has proved to be clearly 
superior to the simplest of all, the 'friend-of-friend' (fof) 
algorithm ( Dav is et al. 19 85 ) . This algorithm links two par- 
ticles in the same group if their distance is less than a certain 
linking length. FOF is very simple, and has the advantage of 
depending on only one parameter, the linking length. This 
is usually expressed as a fraction of the mean inter-particle 
separation. Groups of particles are thus identified which are 
bound by a density contrast 
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For the case where the agglomeration is well-modelled by 
an isothermal sphere, the average density is three times the 
threshold density, so picking a link-length parameter b = 0.2 
identifies groups that are at average densities approximately 
180 times the mean density ( |Cole fe Lacey 1996| l. So, as a 
first step, we detect halos in every output timestep of the 
simulation with this algorithm and create a linked list of par- 
ticles in each detected halo. In practice, we allow 6 to evolve 
with time such that the FOF picks up structures at overden- 
sity around 200 times the critical density. It has been shown 



(KCDW) that the minimum size of groups that are actu- 
ally dynamically stable in numerical simulations is around 
ten particles. We employ, in our fiducial model, a more con- 
servative minimum particle number of twenty. This results 
in FOF groups of minimum mass 1.65 10 n MQ. We find a 
total of 23 000 such groups in our final simulation output 
'snapshot', of which 3000 are not dynamically stable (see 
section l2.5t . 



2.3 Properties 

For each halo, we compute and record: 

• the virial mass of the halo, M200 which is assumed to 
be equal to the group mass, i.e. the number of particles 
multiplied by the particle mass. 

• the position of the centre of mass of the halo, and its 
net velocity. 

• a, b, c, the principal axis lengths of the mass distribu- 
tion. 

• the thermal energy of the halo. This is simply the sum 
of the squares of the particle velocities relative to the net 
velocity of the halo, multiplied by the halo mass. 

• the potential energy. This is computationally expensive 
to measure, since it involves a sum over particle pairs. In 
appendix [B] we explain and justify our approximation for 
measuring this quantity. 

• the virial radius -R200, which is deduced from the 
measured mass and potential energy, assuming a singular 
isothermal sphere density profile for the halo. 

• the circular velocity of the halo, Vaxi- 



V200 — GM200 1 -R200 • 

• the halo spin parameter, A, defined by 
J\E\V 2 



A = 



GM FOF 



(2.2) 



(2.3) 



where J, _E,and AfpoF are the total angular momentum, en- 
ergy and FOF group mass respectively. 



2.3.1 Halo spin distribution 

In the top left panel of Fig. we show the distribution of 
the spin parameter, for halos in our ACDM, A'part = 256 3 
simulation. We also plot the empirical lognormal distribu- 
tion used by |Mo, Mao, fc White (1998| , which is found to 
be a good approximation to the distributions obtained from 
several other investigations (eg. IBarnes fc Kfstathiou~ 1987 
|Cole fc Lacey 19960 . Our distribution is slightly skewed rel- 
ative to this model, with a tail down to low spin parameter, 
but it is peaked at A m 0.04, close to the 'fiducial' median 
value of 0.05. Our median value of A is 0.038 for bound halos 
(see section \2. 51 . and the dispersion in In A is 0.52. This is 
similar to results in these other investigations. 



2.3.2 Halo mass function 

In the top right panel of Fig. we show the halo mass func- 
tion, ie. the number of halos as a function of virial mass. We 
compare this measurement with two analytic predictions. 
The Press-Schechter approach (IPress fc SchechterT974l 
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Figure 1. Distributions of several halo properties. In each plot, the solid black line is for all groups identified by the friends-of- friends 
method, the dotted line is for halos with negative binding energy (ie. bound objects), and the dashed line is for unbound systems. The 
top left panel shows the dimensionless spin parameter, compared to the empirical (log-normal) function of |Mo, Mao, fc White (1998} 
(dashed line, with average A = 0.04, with arbitrary normalization). On the top right is the mass function (ie. distribution of halo virial 
masses), compared to Press— Schechter and Peaks model predictions (short- and long-dashed lines). Bottom left shows the formation 
time, defined in 12. 41 Bottom right shows the dependence of formation time on the halo mass, showing the median and ±1-<t spread of 
the distribution. 



uses a model for spherical collapse of initial perturbations 
to predict the final distribution of halo masses in a given 
cosmology. This is plotted as a short-dashed line in Fig. 



The peaks formalism of |Bardeen et al. ( 1986| > uses the 
general properties of Gaussian random fields, assuming ini- 
tial peaks will be associated with future halos. Predictions 
for the number density of peaks, as a function of peak height 
l |Lacey fc Silk 1991| IGuiderdoni et al. 1998L thus translate 
to a number density of halos as a function of virial mass. 
The prediction from this model is shown by the long-dashed 
line in Fig.0 



It can be seen that neither these two models provides an 
accurate fit over the whole range of halo masses, with Press- 
Schechter overpredicting the number of massive halos, and 
the peaks formalism overpredicting the mass function at the 
light end, below 10 12 solar masses. 



2.4 Building the merger tree 

In hierarchical models of structure formation, small objects 
form first, then merge together into more massive ones, or 
accrete matter from the background. These processes can be 
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Figure 2. An example merging history for a dark matter halo with mass 1.9 10 14 M Q at rcdshift zero. We show all the progenitors of 
the halo as circles with radii representing their virial radii on the same scale as the axes. 



represented with a tree, the merging tree of the dark matter 
halos. Our aim is to follow the evolution of halos detected 
with the group finder in the simulation outputs, and to save 
in a tree structure all the processes of accretion and merging 
between those objects, before subsequently modelling semi- 
analytically the evolution of baryons within these halos. This 
hybrid approach has the advantage that it allow us to use 
all the spatial information from the simulations as inputs to 
our models. 

Considering a halo hi at one timestep, we will identify 
a halo h,2 at the previous timestep as a progenitor if at least 
one particle is common to both halos. 

Having compiled a progenitor list for each halo, we com- 
pute the mass fraction of the halo coming from each pro- 
genitor. This gives the basic structure of the tree. This is a 
non-binary tree, in the sense that each halo can have several 



'sons'. In this general case, the sum of the masses of a halo's 
progenitors can be greater than the mass of the halo itself 
(fragmentation, evaporation, tidal stripping of a progenitor) 
and differs from the implementation of KCDW as each halo 
can have more than one descendent. 

It should be noted that an alternative method has been 
employed by other groups (SP99; CLBF). They have cre- 
ated a merger tree using random realizations of the Press- 
Schechter formalism. Their approach has theoretically infi- 
nite resolution, as the merging history can always be fol- 
lowed back to arbitrarily small progenitors. Our method, 
although resolution limited, has the advantage that we have 
full spatial information about the halos and their progeni- 
tors, and we can look at the 'genuine' history of a halo rather 
than a statistical realization of that history. 

Having defined a list of progenitors for each halo, we 
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can start to examine the history of how each halo acquired 
its present day mass. One particularly simple statistic that 
has been employed to parametrize this history is the forma- 
tion redshift, Z{ of the halo, defined as the time at which 
one half the halo's mass was present in a single object 
l |Lacey fc Cole 1 993 1 . In practice this is calculated by re- 
cursively following the largest progenitor of the halo back 
through the tree, until it has one half the present-day mass. 
If a halo has formed by fragmentation of another halo of 
mass greater than the current mass, we equate the frag- 
mentation time with the formation time. In the bottom left 
hand panel of Fig. Q we plot the distribution of formation 
times, for halos in the final (z = 0) output of our ACDM 
simulation. 

In the bottom right panel of this figure, we show the 
median and scatter (±l-a) of the formation time as a func- 
tion of halo mass, for only those halos with negative binding 
energy. There is a clear trend that more massive halos tend 
to have formed more recently, as expected given the results 
of |Lacey fc Cole ( 1993| > and others. This relationship breaks 
down below « 3 10 M©, showing that, although we resolve 
halos of this size, we cannot resolve their progenitors with 
half this mass, so the formation time becomes the time when 
the halo crossed the mass threshold to become detectable. 

In Fig. [5] we present the merging history of a typical 
massive halo in the final simulation output. We show the lo- 
cations and sizes of all the halo's progenitors in the preced- 
ing timesteps. This clearly shows how virialized halos start 
to form in filaments, and then fall down these structures 
towards the centre, feeding the central massive halo. 

2.5 False halo identification 

No group finder is totally accurate, as it is always possible 
to miss groups or to identify systems that are not in fact 
bound. Having determined the energy of our halos, we find 
that many, especially the lighter ones, have positive total 
energy, ie. they are not bound systems. Leaving these ha- 
los in the sample can seriously skew the spin distribution, 
away from its usual log-normal behaviour. This effect is il- 
lustrated in the top left panel of Fig. where we plot the 
distribution for all the halos, and separately for those with 
positive and negative binding energies. It can be seen that 
excluding halos with positive energy immediately leads to 
an huge improvement in the symmetry of the distribution 
for high spins. 

From the other panels of this figure, we can see that the 
mass function of these halos is steeper than that of 'good' 
halos, and that the majority of them have extremely re- 
cent formation time. The recent formation time in particu- 
lar suggests that the majority of halos with positive binding 
energies are transient collections of particles that have been 
misidentified as bound systems by the group finder. This 
could possibly be remedied by using a different group-finder 
(for instance SO, the spherical overdensity algorithm), or by 
rejecting the 'unbound' particles which contribute most to 
the halo thermal energy. 

However, it must be remembered that some of these 
halos could be made from pairs of similarly-sized objects 
that are currently going through a merger event and have 
not yet relaxed to virialization. In this case, we expect the 
kinetic energy to be rather higher than the virial theo- 



rem would predict, and the halo may have positive bind- 
ing energy. We would not want to reject these halos. Those 
objects with a longer formation time may be smaller ha- 
los in high-density regions that appear disrupted due to 
strong tidal forces, probably shortly before they merge with 
a more massive companion. Clustering studies of these ob- 
jects (Blaizot ct al. 2002) show that they are highly clus- 
tered, and concentrated in the densest regions of the simu- 
lation. 

2.6 Baryon tree 

Although we preserve all the information about particles 
shared, and hence can reconstruct the halo merging tree 
in full detail, for practical purposes we will be interested 
in following only the bulk of the dark matter. Two-body 
(resolution) effects on small scales imply that the exchange 
of one or two particles is not a guarantee that two halos 
are physically related - in a higher resolution simulation, 
this exchange might not occur. For tracing the baryons and 
galaxies between timesteps, we will use a restricted version 
of the tree where we only consider merging between halos 
(i.e. objects with more than 20 particles). Ideally, one would 
like all the mergers to be binary (as in Monte-Carlo gener- 
ated merging trees) , because the order in which halos merge 
might slightly change the properties of the final halo. In 
practice, due to finite time resolution, it is not possible to 
avoid multiple mergers between two outputs of the simula- 
tion, especially for the most massive halos at low redshifts. 
However, we find that most of our halos do not undergo a 
merger between timesteps, and of the merger events that do 
take place, those where three or more progenitors can be 
identified in the preceding timestep are less numerous by a 
factor ten than the number of binary mergers (N mm l999ll . 
We will return to a more detailed study of time resolution 
effects in section |9~^1 but in the mean time we assume, for 
simplicity, that halo mergers always take place exactly in 
the middle of two timesteps. 

When considering a halo, we sort the list of all its pro- 
genitors in decreasing order of contribution to its mass. For 
each of these halos, we then see whether our halo is the 
main (ie. most massive) son. If it is, it is assumed that the 
whole baryonic component (discussed in the next section) 
of the father is transferred to the son. In this simplified de- 
scription, then, each halo has zero (it just formed) or several 
progenitors, and zero or one descendent. We note that this 
simple tree is very similar to KCDW's as their halos are also 
allowed to have many progenitors but only one descendent. 



3 BARYON COOLING IN HALOS 

Having compiled a list of halo properties based on their dark 
matter content, and constructed a merger tree from the ex- 
change of Af-body particles, we come to the point where 
we use semi-analytic recipes to add a baryon component to 
these dark matter halos. The key assumption is that the 
overall behaviour of the matter distribution is determined 
by the dark matter density field, and that the baryon field 
thus represents only a small perturbation to the total den- 
sity. Secondly, it is assumed that gravity is the only force 
operating between halos, so at a given timestep, each halo 
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algorithm, or extreme tidal disruption. We do not allow cool- 
ing in these halos, or in halos with spin parameter A > 0.5, 
the threshold for complete rotational support. 

The characteristic cooling time for the hot g £LS , ELS Et 
function of radius, is given by 

3 nm p kT(r) 

c( } ~2p hot (r)A(T)' {6 - 2) 

where /J,m p is the molecular mass of the gas, phot the hot 
gas mass density, and A(T) is the cooling function, which 
depends on the metallicity and temperature of the gas. We 
use the model of [Sutherland fc Dopita (1993} to find A for 
our halos. We will make the assumption that the gas has 
the same temperature throughout, which is known to be a 
good approximation for clusters, at least outside the central 
cooling flow region (eg. |D avid, Jones, fc Forman 1996) , This 
isothermal approximation, applied to the equation of hydro- 
static equilibrium, results in halo temperature (in Kelvin); 



2k 



(e^t) 2 - 



This latter result comes from fixing [i by assuming the gas is 
totally ionized, consisting of one part helium to three parts 
hydrogen (by mass). We will also assume that the dark mat- 
ter distribution follows an isothermal profile, in which case 
V c , the halo circular velocity, is identical to V200, defined in 
section 12.31 

Given the cooling rate, we can calculate the mass of gas 
cooled during a time interval At: 

m coo i = / p ho t (r)d 3 r (3.4) 
Jo 

In eauation l3.4l r max is the minimum between the infall ra- 
dius ri n f, the cooling radius r coo i, (defined respectively as 
the radii beyond which the gas has not had enough time 
to fall/cool onto the central galaxy in At), and the virial 
radius i?200- We now need a model for the distribution of 
the hot gas in the halo. Observationally, measurements of 
profiles of X-ray gas are often translated into a /3-model 
ICavaliere fc Fu sco-Fcmiano 1976), where the density fol- 
lows 

/ 2 S -3/3/2 

Phot(t) =pa{t) f 1 + ^2 J ■ (3.5) 
Several numerical studies (Eke, Navarro, & Frcnk 1998 



identified in the dark matter field can be considered com- 
pletely independent from all other halos in the simulation. 

3.1 Keeping track of baryons 

When a halo is first identified in the A-body simulation, we 
assign it a mass of hot gas consistent with the primordial 
baryon fraction of the Universe, ie. 

Afhot=M 200 ^ (3.1) 

Q,b is a free parameter in our model, as the power spec- 
trum used to construct the simulations has no dependence 
on the baryon fraction. For the results presented in this pa- 
per, we adopt a value of Sis = 0.02ft -2 , suggested by ob- 
servations of the deuterium abundance in QSO absorption 
lines ffytler, Fan, fc Buries 19961 - 

As the halo subsequently accretes mass, we increase its 
stock of hot gas, assuming that the new mass accreted has 
the same primordial baryon fraction. If the halo loses mass, 
for example through evaporation of dark matter due to two- 
body effects, or from fragmentation (when the halo's prin- 
cipal 'son' is less massive than its progenitor), we discard 
some of its hot gas, based on the current baryon fraction of 
the hot gas in the halo. 

We will see in section 14.21 that galaxy formation pro- 
cesses can also lead to the ejection of gas and metals from 
the halo. Whilst the physical processes are difficult to model, 
it seems reasonable to assume that a certain portion of this 
material will remain in the local environment of the halo, 
and may later be re-accreted. For this reason, we keep track 
of the mass of gas, metals and dark matter that have been 
lost from the halo through these processes, storing them in a 
halo 'reservoir'. It will be assumed that, when the halo subse- 
quently accretes dark matter from the background, some of 
this dark matter will be 'primordial' (ie. with baryon frac- 
tion Qb, and zero metallicity), and the rest will have the 
same baryon fraction and metal content as that of the halo 
reservoir, up to the point where the reservoir is fully used up. 
We parametrize this effect with the halo recycling efficiency, 
C, where ( — implies permanent loss of ejected material, 
and £ = 1 is maximally efficient re-accretion, where all the 
reservoir is used up before the accreted dark matter is as- 
sumed to be primordial. 

There have been previous attempts to deal with this 
issue. |Benson et al. (2000} and KCDW both compare the 
standard model of CLBF in which gas is re-accreted after 
the parent halo has doubled in mass, with a model in which 
re-accretion is immediate (enriched gas can cool back onto 
the disk immediately after it has been ejected). Both of these 
models lack strong physical or observational support. SP99 
adopt the extreme case where ejected material is forever lost 
from the halo. This corresponds to our model with £ = 0. 
For the rest of this paper, we will adopt a standard model 
of C = 0.3. 

3.2 Cooling the hot gas 

At each time step we allow the hot intra-halo medium to 
cool onto a disc at the centre of each halo. As described in 
sect ion 1231 there are some halos which have positive binding 
energy, either due to false identification by the group-finding 



|Navarro, Frenk, fc White 1995b} suggest that clusters are 
well fit by a profile with j3 = 2/3, ie. a non-singular isother- 
mal sphere with core radius r c . Here, we assume that the 
hot gas is in hydrostatic equilibrium within the virial radius 
of our isothermal dark matter haloes, which yields a hot gas 
distribution very similar to that of the /3-model. The cen- 
tral density po of the gas is computed at each timestep by 
integrating over the density profile given by the hydrostatic 
solution and equating the result to the current mass of hot 
gas in the halo (it is therefore time dependent). We can also 
define a 'time integrated cooling radius', r^ ol , as the radius 
that would contain the current cold mass of gas, if all the 
gas (hot and cold) were distributed according to this profile. 
This radius is close to the one used in pure semi-analytic 
models to compute the quantity of gas which has cooled 
since a halo has formed. 



GALICS I 9 



3.3 Overcooling 

A key problem with semi-analytic models is the overpredic- 
tion of the bright end of the galaxy luminosity function. In 
the hierarchical framework, this can be interpreted as mean- 
ing that the central galaxies in massive dark matter halos are 
generally too luminous, either because they have too much 
'fuel' for star formation, or because the merger rate in the 
halo is too high. A number of methods have been employed 
or suggested to remedy this situation: 

• The approach of KCDW has been to prevent gas from 
cooling altogether in halos with circular velocity greater 
than a certain limiting value. Unfortunately, this very sim- 
ple assumption, introduced to account for the observational 
fact that disc galaxies with rotational velocities in excess of 
350kms _1 are extremely rare, has no physical justification. 

• CLBF point out that large halos are formed from small 
halos that have already undergone cooling, and thus they 
will lack the low entropy gas that cools most easily, so one 
should expect the hot gas profile to be rather broader than 
that of small halos, and cooling will be suppressed. Although 
this model has some theoretical basis, there seems to be no 
observational evidence that large clusters have flatter gas 
density profiles. 

• |Nulsen fc Fabian (1995| l point out that if the cooling 
time of some gas is longer than its freefall time to the centre 
of the halo, a cooling flow may occur, in which case the gas 
may remain in a compact cold cloud without forming stars, 
or it may form stars with a non-standard IMF, producing 
an excess of brown dwarfs that will not be visible. 

• |Wu, Fabian, fc Nulsen (2000| l show that the energy in 
the hot X-ray gas is such that there needs to be signifi- 
cant non-gravitational (feedback, AGN) heating of this gas. 
In other words, if we do not include these heat sources we 
may be significantly underestimating the temperature from 
which our hot gas must cool. However, it seems that the 
importance of this effect is likely to diminish with the mass 
of the halo, so this may not in itself solve the problem. 

• Competitive cooling, whereby satellite galaxies can ac- 
crete cold gas as well as the central galaxy in a halo, is gen- 
erally neglected in semi-analytic models, and would again 
help to prevent the assembly of over-luminous galaxies in 
the centres of halos. However, numerical simulations of satel- 
lites moving in a hot intra-cluster medium strongly suggest 
that RAM pressure stripping is the dominant effect (eg. 
IQuilis, Moore, fc Bower (2000| l). 

• SP99 point out that mergers of halos may shock heat 
the gas within, again resulting in the suppression of cooling 
onto the central galaxy. 

The failure to reach a consensus between these dif- 
ferent approaches, and the absence of one well-motivated, 
observationally and theoretically justified model, must re- 
grettably be regarded as one of the failings of the semi- 
analytic treatment, even though it probably finds its roots 
in a poor theoretical understanding of cooling flows. There- 
fore, in what follows, we simply decide to take advantage 
of the observed correlation between AGN and bulge mass 
pa gorrian et al. 1998| l , and prevent gas from cooling in a 
halo when the bulge mass locked in its host galaxies becomes 
equal to 10 11 M . 



We defer a more detailed modelling and investigation 
of the observational consequences of some of the previously 
mentioned models to future work, as this is clearly beyond 
the scope of the present paper. 



4 GALAXIES 

As hot gas cools and falls to the centre of its dark matter 
halo, it settles in a rotationally supported disc. If the spe- 
cific angular momentum of the accreted gas is conserved and 
starts off with the specific angular momentum of the dark 
matter halo l |Mo, Mao, fc White 1998| , we assume it forms 
an exponential disc with scale length ro given by: 

rr> = -^=R 2 oo- (4.1) 
V 2 

As gas builds up over a period of time to form a massive 
central disc we recompute its scale length at each timestep, 
by taking the mass-weighted average gas profile of the disc 
and that of the new matter arriving, whose scale-length is 
also given by 14.11 However, the half-mass radius of the 
disc (see definition below) is never allowed to exceed the 
virial radius of the halo, or to become smaller than the half- 
mass radius of the bulge, provided a bulge is present. We 
do not recompute the disc scale-length of galaxies which 
become satellites. We also obtain the circular velocity of 
the disc, V c by enforcing rotational equilibrium at the disc 
half-mass radius, including all the mass enclosed within, i. e. 
baryons and dark matter core. We assume that the DM core 
is depleted and/or flattened by a rapidly rotating disc as 
we would otherwise overestimate the amount of dark mat- 
ter present within Milky- Way discs pSimiey fc Evans 2001| > . 
More specifically the fraction of dark matter within the half- 
mass radius of a pure disc /dm is taken to be: 

/DM = 1 -(500TnVi) 2 ' ^ 

so that a typical Milky- Way disc only has 80 % of its DM 
core left. This correction should be considered as a first at- 
tempt to model the complex issue of interaction between 
dark matter and baryons, which we plan to address in more 
detail in future work. Finally, we point out that this deriva- 
tion of disc velocity differs from that of CLBF as these au- 
thors use a NFW profile for the density of the dark matter 
halo and include fully self-consistent disc gravity. 

We will be concerned with a number of properties of 
galaxies, including their sizes, and we will generally use the 
radius containing one half the galaxy mass to define a size, 
since this is more physically relevant than the scale-length. 
For an exponential disc, the half-mass radius is given by 

r 1/2 = 1.68r D . (4.3) 

Galaxies remain pure discs if their disc is globally stable 
(i.e. V c < 0.7 x Vtot where Vtot is the circular velocity of the 
disk-bulge-halo system ; see e.g. [van den Bosch (19981 0, and 
they do not undergo a merger with another galaxy. In the 
case where the latter of these two events occurs, we employ 
a recipe which will be described in section |K| to distribute 
the stars and gas in the galaxy between three components in 
the resulting, post-merger galaxy, that is the disc, the bulge, 
and a starburst. In the case of a disc instability, we simply 
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transfer the mass of gas and stars necessary to make the disc 
stable to the burst component, and compute the properties 
of the bulge/burst in a similar fashion as that described in 
section |K| for galaxy mergers. Bulges are assumed to have a 
density profile given by the |Hernquist ( 1990| > model, 

r \ M r B ia a\ 

^=2^W" (44) 

The half-mass radius is given by 

ri /2 = (l + V2)r B . (4.5) 

The bulges are assumed to be pressure supported with a 
characteristic velocity dispersion a, computed at their half- 
mass radius. The material forming a starburst will be as- 
sumed to concentrate in a dense nucleus at the centre of 
the bulge. We assume that its geometry is also modelled by 
Hernquist's equation, but the scale radius, rs, is given by 



rs 



(4.6) 



where k is a free parameter to which we give the fiducial 
value 0.1. Note that gas does not cool onto bulges or star- 
bursts: after a merger, a new disc is formed around the bulge 
by the cooling gas. 

Given a passively evolving galaxy, there are four pro- 
cesses that determine its make-up in terms of gas, stars, 
and metals, these being cooling, star formation, feedback 
and metallicity. Although the metallicity of a galaxy varies 
more slowly, as it takes time for stars to process hydrogen 
into metals and release it back into the gas content, the other 
three can vary on quite short timescales, and are highly cou- 
pled. For this reason, we use an adaptive time resolution for 
calculating all quantities, breaking each simulation timestep 
(or each interval between mergers, whichever is shorter) into 
multiple substeps. Essentially, this fine resolution is a brute- 
force numerical way of solving the differential equations that 
relate cooling, star formation and feedback. 

We now describe the interlinked processes responsible 
for making galaxies. 

4.1 Star formation 

At each timestep, a certain amount of the cold gas in a 
galaxy is allowed to form stars. The rate of this star forma- 
tion is given by 

M co i d 



** = 



(4.7) 



(3 tdyn 

tdyn is the dynamical timescale of the galaxy, for which we 
use the time taken for material at the half-mass radius to 
reach either the opposite side of the galaxy or its center for 
a disc and a bulge respectively, and is given by: 



idyn = ri/2 



for discs 
for bulges 



(4.8) 



where V c is the circular velocity of the disc, and the ma- 
terial is assumed to have purely circular orbits, and a is 
the bulge velocity dispersion, where we assume the matter 
in the bulge has only radial orbits. f3~ 1 is the star forma- 
tion efficiency. Based on data for a sample of bright galaxies 
( Kcn mcutt, Tamblyn, fc Uongdon 1994) it has been shown 
by |Guiderdoni et al. (1998} that a value /3 w 50 is capable 
of reproducing the mean Roberts times jRpbcrts 1963J for 



discs. We therefore use (3 = 50 as the fiducial value of this 
parameter. One can wonder why we define a star formation 
rate for the bulge as gas cannot cool on this component and 
stars form in the burst component during a merger. This is 
simply because stars eject gas and metals back in the ISM 
which can be trapped in the bulge triggering star formation 
there. 

As mentioned earlier, starbursts have the same geom- 
etry as bulges, so their star formation law is the same. 
The only difference is that the radius of the burst being 
smaller, the characteristic timescale for star formation is 
much shorter, giving birth to an almost instantaneous 'burst' 
of star formation. 



4.2 Feedback 

For every mass of stars formed, a certain fraction will be 
contained in massive, fast-evolving stars that quickly be- 
come supernovae. When a supernova occurs, some of the 
energy released will be imparted to the ISM, 'blowing off' 
some of this gas. If this ejected gas has high enough kinetic 
energy, it may even be blown out of the host halo altogether. 
This process thus inhibits future star formation, and hence 
is termed feedback. 

The feedback is given by (ISilk 2 0031: 



2** 



(4.9) 



where e w is the efficiency of the supernova-triggered wind 
which is proportional to v^ ac and depends both on the 
porosity of the ISM (see |Silk (2001) for details) and 
the mass-loading factor. This latter accounts for entrain- 
ment of interstellar gas by the wind and can be con- 
sidered as a free parameter whose value is around 10 
l |Martin, Kobulnicky, fc Heckman 2002) . This is modelled 
by our efficiency parameter e, and as a result, the outflow 
rate given by equation 14.91 for a starburst is of the order of 
the star formation rate. Note that in the previous equation, 
rjsN is the number of supernovae per unit star-forming mass, 
which is a prediction of the Initial Mass Function (IMF) cho- 
sen, and Esn is the energy of a supernova, assumed to be 
10 51 erg. 

Equation 14.91 is applied to find the fraction of gas in 
the ISM that is lost by the galaxy and ejected in the intra 
halo medium. We then equate the fraction of this gas that is 
completely ejected from the halo, to the galaxy/halo escape 
velocity ratio. The gas ejected from the halo is added to 
the halo reservoir where it may subsequently be accreted, as 
discussed in section |3~T1 

The escape velocity is given by: 



In 2 for discs 

for bulges 
Tn(i?2oo/f) for halos 



(4.10) 



where v is the characteristic velocity (ie. the circular velocity 
for discs and halos, the dispersion velocity for bulges) , and, 
in the case of halo ejection, r is either the galaxy's orbital 
radius in the halo, or, for the central galaxy in the halo, the 
half mass radius of the galaxy. 
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4.3 Metallicity 

The baryonic gas in halos and galaxies, initially composed 
solely of hydrogen and helium, acquires a metal content due 
to the processing of these light elements by the stellar popu- 
lation, and the subsequent release of the heavy elements syn- 
thesized to the inter-stellar medium at the end of a star's 
life. To track the metallicity, we need two components: a 
model for the rate at which metals are produced inside the 
stars, and a model for the amount of material released by a 
given stellar population. 

Using these models, we are able to use the tree to track 
the metals ejected by stars at each time step, over the whole 
merging history of the galaxy, and let new stars form out 
of the enriched gas, assuming instantaneous mixing in the 
ISM. This is in contrast to the assumption of instantaneous 
recycling, the approach generally used in semi-analytic mod- 
elling (CLBF, KCDW, SP99). Correct modelling of the ejec- 
tion is likely to be especially important for elements such as 
iron, nitrogen and possibly carbon, for which the production 
delay can be significant (eg. 1 Gyr for Fe produced by SN 
la). 

Given a mass M of stars formed at some time to, we 
can calculate the current contribution to the stellar ejection 
during a timestep from t\ to ti as 

j-m(t 2 ) 

AM+ = -M [m - w(m)]4>(m)dm (4.11) 

J m(t!) 

where m(t) is the mass of a star having lifetime t, w(m) 
is the mass of the remnant left after the star has died, 
and <f>(m) is the IMF. Our fiducial model uses a Kennicutt 
IMF, and in Paper II we will explore the effects of using a 
Scalo or Salpeter IMF (see IKennicutt "19831 |Salpeter 1955| 
IScalo 1986H . For all IMFs considered, we form no stars 
lighter than O.IMq or heavier than 120Mq. 

This formula can be extended to include the entire star 
formation history, producing a rate of increase in the gas 
mass due to ejection from the stellar population, 

E(t) = / **(t-t m )[m- w(mj\<t>{m)dm (4.12) 

J m(t) 

where t m is the lifetime of a star of mass m (ie. the inverse 
of m(t)). 

Equation l4.12l can be adapted to predict the amount of 
metals produced, by making the replacement: 

[m — w(m)] — > [m — w(m)]Z co id (t — t m ) + mY z (m) (4.13) 

where the first term on the right hand side represents the re- 
introduction of the metals that were originally in the stars 
when they formed, and Yz(m) is the fraction of the ini- 
tial stellar mass transformed via stellar nucleosynthesis into 
metals, known as the stellar yield. Our models for the rem- 
nant mass and the yield come from the STARDUST code of 
DGS, and are calculated self-consistently with the models 
for spectral evolution we will discuss in section |S] 

Throughout this work, we assume chemical homogene- 
ity (instantaneous mixing), such that outflows caused by 
feedback processes are assumed to have the same metallicity 
as the inter-stellar medium, though in reality the material 
in the outflow could be metal-enhanced < |Pagel 1998^ . 



5 MERGING 

In the hierarchical picture of structure formation, mergers 
are clearly crucial in understanding the properties of galax- 
ies observed at the present day. In our models, mergers are 
responsible for the formation of massive spheroids, and are 
assumed to trigger starburst activity. The detailed modelling 
of the physics of individual mergers between galaxies is cur- 
rently beyond the scope of what can be achieved in a cos- 
mological simulation. Only the highest resolution A/-body 
simulations of |Springel et al. (2001} are able to dynamically 
follow the merger of 0.2 L* dark matter sub-halos within 
a single cluster. Therefore, the best we can hope for is to 
approximately model the rate and effects of merging in a 
global sense. To deal with the rate of merging, we model 
two effects, the gradual tendency of satellite galaxies to lose 
orbital energy and sink towards the centre of the cluster 
potential well, and the likelihood of occasional encounters 
between satellite galaxies in a cluster, based on probabilis- 
tic, cross-section arguments. 



5.1 Halo mergers 

We identify mergers between dark matter halos using the 
halo tree, as described in section l2~4l When a merger occurs, 
the properties of the dark matter itself are obtained directly 
from the properties of the new halo in the ./V-body simula- 
tion. We apply recipes, however, to deal with the baryonic 
component. Firstly, the ICM, ie. the hot cluster gas, of the 
progenitors are added together proportionally to the dark 
matter mass which ends up in the descendent halo (the con- 
stant of proportionality being the current baryon fraction 
of each progenitor), and given the virial temperature of the 
new halo. Then, the properties of the two halo reservoirs 
are simply added. Note that the new halo contains all the 
galaxies that were present in its progenitors (even if frag- 
mentation occurs and the halo is less massive than all/some 
of its progenitors) and that the fraction of the progenitors' 
ICM which does not end up in the new halo is put in its 
reservoir, ensuring conservation of metals. 

When two halos merge (sometime in between two out- 
put times) , there is a discontinuity between the centre of the 
new halo and the centres of mass of the progenitors (linearly 
extrapolated from their positions and velocities in the pre- 
vious timestep). We measure this 'jumping' distance, Rj, for 
each of the progenitors directly from the iV-body simulations 
and use it to assign each galaxy an orbital radius in the new 
halo, using the cosine law: 

r nC w = yj rl ld + Rf - 2r oM Rj cos 6, (5.1) 

where r id is the orbital radius of the galaxy in the progeni- 
tor halo and r ne w is its orbital radius in the new halo. This 
assumes that the galaxy was initially at an angle 6 to the 
vector joining the centres of the two halos. We select cos 6 
randomly from a flat probability distribution between — 1 
and 1, taking account of the spherical symmetry. 

To illustrate how this works in practice, we consider two 
cases. The first is a merger between a small halo and a much 
larger halo. It is clear that the centre of mass of the large 
halo is perturbed only slightly by the encounter. Rj is thus 
small, and the galaxies that were previously in the large halo 
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have r ne 



Told 



For the small halo, R\ is close to the virial 



radius of the large halo. Thus r ncw ~ Rj for the galaxies 
that were in this halo, in other words they are placed close 
to the virial radius of the new halo. For a collision between 
equal mass halos, Rj is approximately the progenitor virial 
radius, so the central galaxies are placed at this distance, 
while non-central ones are placed randomly throughout the 
new halo. 

Any galaxy whose orbital radius after the merger is 
less than its own half-mass radius becomes the new cen- 
tral galaxy of the halo. If there is more than one of these 
objects, they are merged with each other in order of de- 
scending mass. It may well happen that there is no central 
galaxy after a merger, in which case we cool gas onto the 
galaxy closest to the halo centre. Although this scheme is 
really only a naive geometrical model, we use it to approx- 
imately replicate the sort of scatter in galaxy positions we 
expect to see when halos collide. It has the advantage that 
it reproduces the desired behaviour for the extreme cases of 
equal or very unequal mass mergers in a natural, continuous 
way without free parameters. 



5.2 Dynamical friction 

Dynamical friction causes satellite galaxies to sink gradu- 
ally to the centre of their host halo, resulting in a merger 
if there is already a galaxy at the centre. After a merger 
between halos, the galaxies are given orbits as prescribed in 
section |5~T1 In subsequent timesteps, the radii of the orbits 
are decreased by an amount obtained from the differential 
equation (Binn ey" fc Tremaine 1987| l: 

* _ .428%2MnA (5.2) 
dt V c y ' 



where r is the orbital radius, m sa t the mass of the satellite 
(including its tidally stripped dark matter core, see next 
section), V c the circular velocity of the halo and In A the 
Coulomb logarithm, approximated by 



A = 1 + 



V m sat / 



(5.3) 



(see eg. Mamon 1995J. Calculating the amount of orbital 
decay at each timestep seems far more natural than the 
approach favoured by other authors (CLBF, KCDW) who 
assign each galaxy a dynamical friction timescale, and as- 
sume the merger occurs after that fixed time, unless a major 
merger occurs, at which point this timescale is recomputed. 
By modelling the true radial coordinate, albeit in a naive 
way, we automatically take into account possible evolution 
of the halo properties, and avoid having to differentiate be- 
tween 'major' and 'minor' mergers. 

As soon as the orbital distance of a galaxy becomes 
lower than the sum of its half-mass radius and the half- 
mass radius of the central galaxy in the halo, it is assumed 
to merge with this central galaxy. 

Other workers (CLBF, KCDW, SP99) have taken into 
account the fact that galaxy orbits are not in general cir- 
cular. This fact alters the extent of dynamical friction, and 
|Navarro, Frenk, fc White ( 1995a} have shown that applying 
the formula for circular orbits results an average overpredic- 
tion by a factor two of the time taken for an infalling galaxy 
to reach the centre. A proper consideration of this effect re- 



quires a modified form of the dynamical friction as a function 
of orbit eccentricity l |Lacey fc Cole 1993} , and an assump- 
tion for the distribution of ellipticities. However, careful 
examination of |JNavarro, Frenk, fc White 199 5a (especially 
their figure 8), shows that one can obtain the true dynami- 
cal friction time more accurately by simply halving the time 
for circular orbits, than by using the actual orbital eccentric- 
ity and applying the modified form of |Lacey fc Cole 1993| 
We thus take into account the effect of non-circular orbits by 
multiplying the right hand side of equation 15.21 by a factor 
two. 



5.3 Tidal stripping 

When a merger between halos occurs, the new satellite 
galaxies will, to some extent, retain the dark-matter cores 
of their previous host halos, and this dark matter will con- 
tinue to dominate their dynamics. However, as the galaxy 
slowly descends to the centre of the new host this sub-halo 
will be stripped by tidal interactions. Eventually, as the 
mass of the remnant becomes comparable to the baryonic 
mass of the galaxy itself, it will no-longer dominate over the 
effects of the self-gravity of the galaxy, and the dynamics 
may undergo a serious change. Following SP99, we assume 
that the sub-halo is stripped down to a radius rt such that 
pc(rt) = Ph{ro), where p c is the core density and Ph(ro) is 
the halo density at the orbital radius of the remnant (for 
isothermal spheres this density is simply proportional to the 
mean density inside the same radius). Note that we do not 
attempt to model the stripping of baryons by tidal effects, 
only that of the dark matter. 



5.4 Satellite-satellite mergers 

Satellite mergers are assumed to occur through direct colli- 
sions between galaxies in clusters. In the absence of detailed 
numerical simulation of the precise positions and velocities 
of the galaxies within a halo, it is appropriate that some sort 
of cross-section argument be used to determine the merging 
rate. There is no detailed study of how such cross-sections 
behave in a realistic halo environment; a first analytical 
es timate by |Mamon (1992^ has been followed numerically 
by |Makino fc Hut (1997) , with |Mamon (2000} showing that 
these results are in good agreement. We assume the proba- 
bility of a galaxy having a merger in time At is given by: 



At 



(5.4) 



where r is the merger timescale, whose dependence 
on galaxy and halo characteristics is parametrized by 
Maki no fc Hutl as: 



«gal_ 

ri/2 



(5.5) 



In this formula, iV g is the number of galaxies in the halo, 
and Ugai is taken to be a mass-weighted average of the disc 
rotation speed and the bulge dispersion velocity. A direct 
transformation of the results of |Makino fc Hut (1997) leads 
to a value of ip = 0.017 sGyr _1 km~ Mpc, which is the value 
we shall take for our fiducial model. However, their work as- 
sumes all galaxies are spheroids of the same size and mass, 
which is not the case in our models, and in Paper II we 
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Figure 3. Our algorithm for calculating the fraction of disc ma- 
terial remaining in the disc after a merger, as a function of the 
mass ratio of the two progenitors. 



will investigate the effect of altering ip from this standard 
value. We are aware that our extrapolation of these authors' 
results is quite crude and should be, in the best of cases 
taken as a rough estimate of the importance of satellite- 
satellite mergers. Finally we note that the results obtained 
by |Springel et al. (2001} for sub-halo merging in dark mat- 
ter clusters support the view that this kind of mergers is 
negligible in such environments. 



5.5 Post-merger morphology 

In our model, mergers and disc instabilities are both respon- 
sible for creating galaxy bulges, and as such are drivers of 
morphological evolution. In the case of a disc instability, 
the mass transfer between galaxy components has already 
been discussed in section 01 and therefore we will only be 
concerned by mergers in this section. The standard method 
for modelling a merger-driven morphology (SP99, CLBF, 
KCDW) has been to take the ratio of the masses of the 
two objects merging, and add the stars of the lighter galaxy 
to the disc of the heavier one if the mass ratio is less than 
/bulge ~ 0.3, or destroy the disc and form a bulge if the ratio 
is higher. 

It has been shown jWalker, Mihos, fc Hernquist 1996| > 
that a galactic disc can be completely disrupted by an en- 
counter even when the interloper is rather less massive than 
the disc itself, and this simple method reproduces this be- 
haviour. However, it does not allow for any intermediate 
behaviour - there is a sharp cut-off at /bulge between totally 
disrupting, or not disrupting, the morphology of the galaxy, 
which seems quite false. 

Our model is constructed in terms of the function X. In 
general, this can be any function of the ratio of progenitor 
masses that maps the range one to infinity onto the space 
to 1. The step function used be other workers is one example 



of such a function, but we will consider a smoother function, 
drawing inspiration from Fermi-Dirac statistics, and defining 



X(R) 



(5.6) 



This function is shown in Fig.|^] for the value \ = 3.333. 
In equation 15.61 R represents the mass ratio of the heavier 
to the lighter progenitor, and \ i s the critical value of that 
ratio, ie. the value that gives X = 0.5. Since step functions 
with /bulge ~ 0.3 have been found to give good results in the 
past, we set \ — 1/0.3 as our fiducial value. In Paper II we 
will explore the effects of varying \- 

We then define two matrices, 



A s tar — 



X 

1 

1- X 1 

X 



l-X 1 1 



(5.7) 



(5.8) 



Considering the vector V = (D, B, S) where D, B, S refer 
to the mass in the disc, bulge and starburst respectively, we 
have: 



V 11 



A V 1 

"gas v ga 



+ A V 

I "gas v gas 



(5.9) 



and similarly for the stars. 

In other words, during a merger a fraction X of gas and 
stars originally sitting in the disc remain in the disc. The 
rest of the gas from the disc goes into the starburst, as well 
as its stars. Any stars that were already in the bulge stay 
in the bulge, but all its gas is put in the starburst. All the 
material (gas and stars) that was originally in a starburst 
remains in that starburst. Note that gas is never added to 
the bulge in this process, but that a small amount of gas is 
generally present in the bulges, coming from the evolution 
of its stellar content. 



5.6 Properties of merger remnant 

Having ascertained how much mass goes where, we must 
decide on the properties of the resulting galaxy after the 
collision, ie. its rotation/dispersion speed and size. For this, 
we adopt a similar model to that of CLBF. 

Under the virial theorem, the total internal energy is 
given by E = —T where T is the internal kinetic energy. 
Applying conservation of energy, 



Ti + T 2 - E- a 



(5.10) 



where Eint is the interaction energy of the collision. For a dy- 
namical friction encounter as well for satellite-satellite merg- 
ers, we use the total orbital energy, 

Gm\rri2 



ri + r 2 



(5.11) 



(ri and r 2 are the half- mass radii of the two progenitors). 
Note that mi and m 2 stand for all the mass (dark matter 
core and baryons) contained in each progenitor respectively. 

These energy considerations, coupled to mass conserva- 
tion enable us to obtain r ncw , the half mass radius of the 
merger product: 
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(mi + m 2 ) 2 

-^mt I 1. I 2 

" 0.4G "T" r x T r 2 



(5.12) 



where the factor 0.4 depends (not sensitively) on the exact 
density profile (Hernquist in our case). As for the pure disc 
case, we do not allow r llcw to become bigger than 7?200- We 
then assume that r nC w is the half mass radius of the bulge 
and hypothesize that the largest disc is the most likely to 
survive the merger. Following the method presented in sec- 
tion 01 we can calculate the characteristic velocities of our 
galaxy components. In particular, we use equation 14.61 to 
calculate the radius of the starburst after merging and as- 
sume the starburst has the same velocity dispersion as the 
bulge. 

In the case of a satellite-satellite merger, the orbital 
radius of the merger remnant is obtained through a simple 
mass-weighting of the orbital radii of the two progenitors. 



6 GALAXY LUMINOSITIES 

Having decided where the baryons are, and what their phys- 
ical properties are, we apply a set of models to calculate 
the amount of light they produce. Luminosities in visual 
and near-infrared (rest-frame) wavelengths are calculated 
from stellar synthesis models, which take into account the 
metallicity and age of the stellar population. We then add 
geometry- and metallicity-dependent models for the absorp- 
tion and re-emission of this light by the dust and gas in the 
inter-stellar medium. 



6.1 Stellar Spectra 

For each time output of the simulation, we compute the 
stellar contribution F£(t) to the galactic flux at time t, which 
can be written: 



9i.it -T)4>(m)fx(jn,T, Z^dmdr , (6.1) 



where f\(m,r,Z+) is the flux at wavelength A of a star 
with initial mass m, initial metallicity Z*, and age r (ie. 
t — corresponds to the zero age main sequence, and 
fx(m, t, Z+) = if r > t(m)). <f)(m) is once more the IMF. 
We neglect the nebular component. For f\(m,T, Z+) we use 
the STARDUST model of DGS, and full details can be found 
in that paper. 

We use the full merging history of the galaxy to com- 
pute the stellar spectrum, summing up the contribution to 
the present day spectrum from all the stars formed in the 
previous timesteps in all the progenitors of the galaxy. 



6.2 Dust Absorption 

To estimate the stellar flux absorbed by the interstellar 
medium in a galaxy, one first needs to compute its optical 
depth. As in |Guiderdoni fc Rocca-Volmerange 1987] we as- 
sume that the mean perpendicular optical depth of a galaxy 
at wavelength A is: 



\A V J : 



(Nh) 



2.1 10 21 atoms cur 



(6.2) 



where the mean H column density (accounting for the pres- 
ence of helium) is given by: 



<^H> 



atoms cm 



(6.3) 



lAm p T7(ar 1/2 ) 

and a is calculated such that the column density represents 
the average (mass-weighted) column density of the compo- 
nent, and is 1.68 for discs, 1.02 for bulges and starbursts. The 
extinction curve depends on the gas metallicity Z s accord- 
ing to power-law interpolations based on the Solar neigh- 
bourhood and the Large and Small Magellanic Clouds, with 
s = 1.35 for A < 2000 A and s = 1.6 for A > 2000A (see 
|Guiderdoni fc Rocca-Volmerange 1 987 for details). The ex- 
tinction curve for solar metallicity (A\/Ay)z e is taken from 
|Mathis, Mezger, fc Panagia 1983] 

For the spherical components we use the generaliza- 
tion given by |Lucy et al. (1989^ for the analytic formula 
giving obscuration as a function of optical depth r^ ph 
JOster brock 1989 ) to the case where scattering is taken 
into account via the dust albedo, uj\ (we use the model of 
IDraine fc Lee 198 4 for the albedo): 



A A (r) = -2.51og 10 
where 



"A 



1 — u\ + u\a\ 



3 
4t a 



1 



2r x 2 



+ ^2 1 exp(-2r A ; 



(6.4) 



(6.5) 



We assume dust and stars are distributed together uni- 
formely and thus there is no need to average over inclination 
angle, since the bulge is spherically symmetric. 

For discs, the situation is more involved, due 
to inclination effects. We use a uniform slab model 
jGuiderdoni fc Rocca-Volmera nge 1987) for the extinction 
as a function of the inclination angle, 9: 



A A (T,0) = -2.51o gl 



1 — exp(— a\ sec 9) 
a\ sec 9 



(6.6) 



where ax = \/l ~~ ^aT"a- To calculate this quantity, we pick 
a random inclination of the galactic disc to the observer's 
line of sight. For comparisons of the same galaxy at differ- 
ent epochs, and for examining inclination-corrected statis- 
tics (eg. the Tully-Fisher relation), we also compute and 
store the face-on magnitudes of the discs (sec# = 1). 

Our final result is the extinguished stellar spectrum: 



Fx = Ft 



dex[-0.4A A ]. 



(6.7) 



6.3 Dust Emission 



Dust absorption in bulges is spherically symmetric, so the 
emission is trivially calculated, but for discs it is anisotropic, 
so to calculate the total energy budget available for infrared 
emission by dust one must include contributions from all di- 
rections. We use the angle-averaged version of equation 16.61 



Aa(t) = -2.51o gl 



1 - exp(-aA/a:) 
ax/x 



tlx 



(6.8) 



The value of the integral on the right hand side of this equa- 
tion is: 



[l + (a x - 1) exp(-a A ) - a\ Ei(a A )] , 



(6.9) 
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where Ei is the first-order exponential integral. 

Given an unextinguished stellar spectrum F£, the total 
bolometric infrared luminosity of the galaxy is: 



meaning 



fiducial value 



range 



Fx X (1 - dex[-0.4^ A ]) dX. 



(6.10) 



We assume the emission is isotropic, since dust it- 
self is generally optically thin to infrared light (except 
in heavily obscured starbursts). To compute the IR emis- 
sion spectrum, one needs to model both the size distribu- 
tion and the chemical composition of dust grains in the 
ISM. The model we use is based upon the MW model of 
|Desert, Boulanger, fc Puget (1990j l, which includes contri- 
butions from polycyclic aromatic hydrocarbons, very small 
grains and big grains. Our chief refinement to this model 
is to allow a second population of big grains, closer to the 
star forming region, to be in thermal equilibrium at a higher 
temperature for galaxies undergoing massive starbursts. Our 
emission model makes use of the colour/luminosity correla- 
tions observed by IRAS, and is detailed in DGS. We stress 
that a key weakness of this model is that it is based on a 
local sample and therefore assumes that dust properties do 
not evolve with time. 



7 FREE PARAMETERS 

The semi-analytic recipes described here contain a number 
of free parameters that will affect our results. We detail them 
in this section. 

• Qb is the baryon fraction of the Universe. Whilst a high 
baryon fraction can have a serious effect on both the shape 
of the initial power spectrum of density perturbations (eg. 
Sugiyama 1995| and the internal dynamics of halos, provid- 
ing it remains small it can be treated as a free parameter of 
the galaxy formation recipes. Increasing Qb results in more 
fuel for the star formation process, and should thus produce 
brighter galaxies. 

• /3 _1 is the star formation efficiency. This is the param- 
eter that appears in equation 14.71 Increasing f3 will clearly 
result in less cold gas being turned into stars, moving the 
peak in the star formation history closer to the present day. 

• e _1 represents the efficiency of mass-loading during the 
trigerring of a galactic wind by SN explosions (section I4.2H . 
Decreasing e produces more feedback, heating more cold gas, 
ejecting more hot gas from halos, and thus reducing the 
amount of gas that can potentially form stars. 

• ip is the normalization for satellite-satellite merging 
law, equation 15.51 From fMaki no fc Hut 1997t we expect 
its value to be around 0.017, but as stated in section T5.4I 
we will treat it as a free parameter. Satellite-satellite merg- 
ing competes with dynamical friction merging (especially for 
massive galaxies, see IMamon 2000 ) . and energetics of these 
collisions are different, so increasing ip produces more diffuse 
galaxies, slowing the global star formation rate. 

• x parametrizes our merging law in equation 15.61 it is 
the critical mass ratio, ie. the one at which half the disc is 
disrupted. Our fiducial value is 1/0.3, raising \ means that a 
larger disc will be disrupted by a smaller interloper. It there- 
fore affects the morphological mix of our galaxy population. 

• C, is the efficiency with which gas stored in the halo 
reservoir is admitted back to the halo during accretion, as 



n B 


Baryon fraction 


0.02/1" 2 


0.01-0.03 h~ 2 


P 


Inverse star formation efficiency 


50 


10-100 


e 


Inverse of mass-loading for feedback 


0.1 


0.01-1.0 




S-S merging normalization 


0.017 


0.01-0.03 


X 


galaxy merger power-law 


3.333 


2.0-5.0 


C 


recycling efficiency 


0.3 


0.0-1.0 


K 


ratio of burst to bulge radius 


0.1 


0.05-0.2 



Table 2. Summary of free parameters in the models 

explained in section 13.11 Turning £ right down will mean 
that metals and gas ejected from a halo are only returned 
very slowly, restricting the supply of metals to the cooling 
gas. 

• k is used to determine the starburst radius from the 
bulge radius, as shown in equation 14.61 Smaller k produces 
smaller starburst regions, where star formation is therefore 
faster, brighter, and more heavily extinguished. 

These parameters are summarised in table [5] In addi- 
tion, there are various other inputs to the models, including 
cosmology, initial stellar mass function (our fiducial model 
uses the |Kennicutt (1983| l IMF), and dust recipe, that will 
have an impact. We will consider the impact of changing 
some of these inputs, and varying the above free parameters 
within the range quoted in table [5] in paper II. 



8 RESULTS 

In order to test our models against observations, we will use 
a number of constraints in the literature. These fall into two 
main categories; intrinsic properties, ie. relationships that 
apply for individual galaxies such as the Tully-Fisher rela- 
tion or the dependence of metallicity on galaxy mass; and 
global properties, such as luminosity functions and cluster- 
ing statistics. In this paper we will just present work con- 
cerning the local (ie. low-redshift) Universe, and we will not 
use any spatial information. Paper II will include results 
from higher redshift, and the effect on the results of changes 
in the parameters described in section |7| Papers III and IV 
will be concerned with the evolution of many of the statis- 
tics presented here, including the luminosity function and 
predictions of faint galaxy counts. We will also employ the 
spatial information from the dark matter simulations to in- 
vestigate the clustering properties of our galaxies in paper 
V. 

In Fig. 2] we present a slice through the simulation at 
its final output time, with the dark matter density shown 
as a grey-scale. The slice is of side length 100/i _1 Mpc, the 
whole simulation box, but is 10/i~ Mpc in thickness. We 
overplot the locations of bright {Mb < —19.5) galaxies, 
with the sizes of the points representing the galaxy B-band 
magnitudes. One can clearly see the way that the galaxies 
trace the high-density regions of the simulation, with most of 
the light concentrated in filaments and clusters. The power 
of the hybrid approach is immediately apparent when one 
considers that state-of-the-art attempts to model the sub- 
grid physics with smoothed particle hydrodynamics produce 
around 2000 galaxies in a cube with one-third the volume of 
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Figure 4. A slice through our simulation at redshift 0. Dark matter density is represented by a grey-scale, whilst white circles mark 
galaxies, with their sizes depending on the B-band luminosity. Only bright galaxies (Mg < —19.5) are shown. The slice is 100/i -1 Mpc 
on a side (the entire simulation cube) and 10h~ 1 Mpc in thickness. 



ours (Pcarcc ct al. 20011, when we have 30 000 galaxies in 
total in our final timestep, thus giving us access to a much 
broader range of galaxy mass and merging history. 

8.1 Effect of resolution 

The identification of halos with a given cut-off mass pro- 
duces an effective resolution limit on our results. We define 
the galaxy resolution limit as the smallest halo mass times 
the baryon fraction of the Universe. Galaxies of baryonic 
mass less than this exist in our models, but no galaxy of 
mass greater than this could exist in unresolved halos. This 



is thus a 'safe' limit for completeness in our sample, and pro- 
duces a formal completeness limit of galaxy baryonic mass 
m ga i = 2.2 10 10 Af© for the ACDM simulation. 

fn Fig. |S]we show the dependence of galaxy luminosity 
on total baryonic mass for galaxies in our fiducial model, de- 
fined in section |7| The tight relationship between these two 
quantities allows us to estimate the magnitude complete- 
ness limit corresponding to the mass completeness limit we 
just defined. For each waveband, we perform a linear fit of 
magnitude on luminosity, and take the intercept of this fit 
with the line m ga i = 2.2 1O 1O M0. The results of this treat- 
ment for a number of optical and near-infrared filters are 
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Figure 5. Galaxy absolute magnitudes as a function of total 
baryonic mass in the /-band and S-band. The solid vertical line 
shows the effective completeness limit as defined in section 18. fl 
Results for the completeness limit in a number of wavebands are 
presented in table 131 



Waveband 


Completeness limit 


U 


-18.8 mag 


B 


— 18.9 mag 


V 


— 19.7 mag 


R 


—20.4 mag 


1 


-20.9 mag 


.1 


—21.6 mag 


H 


-22.2 mag 


K 


-22.7 mag 


12 (im 


0.41 1O 9 L 


25 fim 


0.25 1O 9 L 


60 (im 


0.79 1O 9 L 


100 /im 


1.34 1O 9 L 



Table 3. Resolution limits in a variety of optical and infrared 
bands (magnitudes for Johnson-Mould filters, \L\ in solar lu- 
minosities for the IRAS bands). Magnitudes are absolute, in the 
Vega system, and calculated using the process described in sec- 
tion |0] 

presented in table El It will be noted that, in the B-band 
the resolution limit is about 1.5 magnitude fainter than the 
value of M t = —20.5 (assuming our value of Hq) from the 
Stromlo-APM survey ( |Loveday et al. 1992) . We are there- 
fore complete to galaxies around one-fourth as luminous as 
an M* galaxy. 

Although we are effectively complete down to these 
magnitude limits, we do not fully resolve the merging history 
of the faintest objects. This is an advantage of the method 
employed by |Benson et al. (2000| l , where halos identified in 
a dark matter simulation have their merging histories com- 
puted with the Press-Schechter formalism with theoretically 
infinite resolution, so even the smallest objects have a fully 
resolved history. The smallest galaxies in our models are in 



halos that have only recently been identified in the simula- 
tions, and hence have a high gas content and star formation 
rate, but low metal and stellar content, and an inevitable 
spiral morphology. We expect the effects of this lack of in- 
formation to propagate over the formal mass completeness 
limit. We consider resolution effects in rather more detail in 
section 



8.2 The halo ICM 

Before examining the galaxies themselves, we will look at the 
overall behaviour of the baryons in the dark matter halos. 
In the left panel of Fig. |SJ we show the ratio of hot gas 
to total baryonic mass for dark matter halos as a function 
of their circular velocity. Note that for isothermal spheres 
there is a one-to-one mapping from circular velocity to halo 
dark- matter mass, V c oc M 1//3 . 
We notice that: 

• the most massive objects have around seventy per cent 
of their baryonic mass in the form of hot intra-cluster 
medium (ICM), with a weak dependence on circular velocity. 

• for objects with lower circular velocity, V c < 
300 kms -1 , the gas fraction starts to fall off rapidly as the 
halo mass decreases, this scale defining the difference be- 
tween a galaxy and a cluster. 

• finally, around 100 kms -1 , the scatter becomes very 
large. This is due to the finite resolution of our approach, in 
that the smallest objects do not have well-resolved merging 
histories, so cooling does not start until the timestep they 
are discovered, resulting in an overestimate of the hot gas 
fraction. The cooling in these objects is rapid enough, how- 
ever, that by the time the mass has increased slightly, they 
have joined the main locus on the plot. 

In the right hand panel of Fig. |S]we show the behaviour 
of the metallicity (defined as the mass in metals relative to 
the mass in gas) of the hot ICM. The resolution effects are 
clearly obvious again at low halo mass, with many objects 
of extremely low metallicity, since their stars are very young 
and have not had time to contribute via nucleosynthesis and 
ejection. Above 100 kms -1 , the metallicity is a decreasing 
function of halo mass, suggesting that larger halos have a 
higher proportion of primordial gas, and galactic winds are 
less important in clusters than in galaxies. This result is 
in contrast to the observational results of |Renzini (1999"j l, 
who finds for the largest clusters (V c > 1000 kms -1 ) that 
the metallicity, as measured from the iron abundance, has 
a constant value of around one third solar. This is between 
five and ten times more metals than the amount we find in 
our halos. 

We will show in section lcT4*l that our galaxy metallicities 
appear quite accurate. Could the missing metals be hiding 
in the galaxies, or possibly in the halo reservoir? We have 
done the same analysis as in Fig.|SJ but added to the ICM all 
the metals contained in the halo reservoir, and all the metals 
locked up in the cold gas inside galaxies. This thus represents 
the maximum metallicity we can attain without somehow 
lowering the mass of non-metals in the ICM. Although this 
has a large effect for small halos, there is little change in the 
metallicity of the largest clusters. We conclude that we are 
not simply 'misplacing' the heavy elements. 



18 Hatton et al. 




There are several potential explanations for this dis- 
crepancy: 



• observations measure only the abundance at the centre 
of the cluster, assuming the hot gas is chemically homoge- 
neous. In fact, metallicity gradients have been observed in 
many clusters (eg. |Koyama, Takano, fc Tawara 1991| for the 
Virgo cluster. IVVhite etHu^T994l for several cooling flow clus- 
ters), and the observations thus yield an upper limit on the 
total metallicity. 

• the metallicity is measured using the iron abundance, 
whereas the metallicity we quote is the total mass of metals 
in the cluster. In fact, our models pay no attention to the 
energy injected to the ISM from type la supernovae (though 
these supernovae do contribute to the metallicity). The over- 
all contribution of these objects is likely to be small, but, 
since these stars contain an iron core, the relative abundance 
of iron can be strongly enhanced, especially from evolved, 
gas-poor, early-type galaxies, such as are found in clusters. 
Thus the simple assumption used to convert Fe abundance 
to total metallicity, Z/Zq = Z Fc /Z^, may be invalid in 
clusters, and will produce an upper limit on the total metal- 
licity. 

• one should also consider the existence of population III 
stars in the early Universe which may have reionized the 
Universe homogeneously and likewise introduced metals into 
the primordial gas. This could be simulated by allowing a 
non-zero primordial metal abundance in our models. 

• the yields we use could be too low especially if the IMF 
is biased towards massive stars in the early stages of cluster 
formation. 

• Although these effects could all play a role in the metal- 
licity underestimate, it seems more likely that the major 
contribution we are unable to model comes from objects be- 
low our resolution limit. Indeed, it is sensible to argue that 
the majority of metals seen in the ICM today came from 
old dwarf galaxies (IRenzini T9991 iRenamTO OO ) . in shallow 
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Figure 7. Mass to light ratios of our halos, in units of the solar 
mass to light ratio. The upper line is for the B-band, the lower for 
the X-band. The lines are the medians, with error bars showing 
the ±1-<t scatter. We only include contributions from galaxies 
with baryonic mass greater than the formal completeness limit 
(section 18.11 . 



potential wells, and we are a long way from achieving the 
resolution needed to see these galaxies. 



8.3 Mass to light ratios 

The mass to light ratio of dark matter halos tells us impor- 
tant information about measuring the dark mass from ob- 
servations, and about which environments are most efficient 
at producing energy. In Fig.|7|we show mass to light ratios, 
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in solar units, for our halos. We only include contributions 
to the luminosity from galaxies above our resolution limit. 
We compare results in two wavebands, the B-band, where 
the energy is dominated by young stars, and the if-band, 
which has been shown to reliably estimate the total stel- 
lar content (Kauff mann fc Cha riot 199§J|. O ne can see that 
the distributions are qualitatively similar in shape and that 
low mass, galactic halos are brighter per unit mass than the 
larger, cluster-sized objects. In the B-band, MW type ha- 
los have a mass-to-light ratio of around 60, rising to 300 for 
the biggest clusters in the sample whereas in the if-bandthe 
mass-to-light ratio is rather constant around 30 for galaxy 
size halos (V c < 300 kms -1 ) and rises to 60 for the biggest 
objects. 

Our failure to fully resolve the merging history of halos 
with circular velocity less than 100 kms -1 makes it difficult 
to draw any firm conclusions about the shape of the mass- 
to-light ratio at the low-mass end. 



8.4 Assorted galaxy properties 

In Fig. |H1 we show scatter plots for some of the intrinsic 
properties of our galaxies as a function of galaxy char- 
acteristic velocity. This is the circular velocity for discs, 
three-dimensional dispersion velocity for bulges, and a mass- 
weighted average of the two for objects of mixed morphology. 

Firstly, we show the total mass of baryons, ie. cold gas 
plus stars, for our galaxies. The main locus of these points is 

5/2 

seen to approximately follow the relationship A/bar oc V GAh , 
although there is a substantial scatter to low baryon mass for 
galaxies with low velocity. This is another example of reso- 
lution effects entering our results, in that these are galaxies 
in small halos that have not yet had time to cool. Com- 
pared to halos, whose virial masses go as the cube of the 
circular velocity, we already see a difference for galaxies: the 
baryonic mass depends more weakly on the velocity. We can 
already see the effects of the cooling and merging, in that the 
galaxy mass and velocity are obviously not scaling linearly 
with those of its parent halo. 

Secondly, we present the metallicity of the cold gas in 
the galaxy. This is a slowly varying function of circular veloc- 
ity. One can contrast this with the results for the halo metal- 
licity presented in Fig|H] although the metallicity of the clus- 
ter ICM decreases with increasing cluster mass, the metal- 
licity of the central object in the cluster increases. There 
is also a clear upper cut-off in metallicity around twice the 
solar value, which corresponds to galaxies having consumed 
almost all their gas and therefore having a very low, if any 
star formation rate at all, and the very few metal rich young 
stars produced fail to enrich the ISM much more. It is worth 
noting that the main locus of our points is around a flat line 
with Z ~ Zq, indicating that for our mass resolution the 
metallicity of the Milky-way is quite typical. 

In the bottom left panel, we show the star formation 
rate as a function of characteristic velocity. We note that the 
majority of galaxies lie close to the line defined by the star 
formation rate being proportional to the velocity squared, 
but that there is a small population with very high star 
formation rates of up to 100 solar masses per year, pre- 
sumably representing galaxies which are undergoing merger- 
driven starbursts. There is also considerable scatter at the 



property 


actual MW 


reference range adopted 


"Igas/mbar 


0.10 


1 ±0.05 


M K 


-23.7 mag 


2 ±0.3 mag 


K 


220 kms -1 


3 ±20 kms" 1 



Table 4. Observed parameters of the Milky Way, used to 
define a MW type galaxy in the simulation. The references 
are: 1. IPrantzos fc Aubert 19951 2. |Kent, Dame, fc Fazio 1991| ; 
3. |Binney &: Tremaine 1987| 



property 


actual MW 


reference 


Simulation 


"Ibar 


0.6 lO n M 


1 


(0.5 ±0.1) 10 n M Q 




2-4M Q yr- 1 


1 


(1.5±O.4)M yr- 1 


'D 


4-5 kpc 


2 


(2.8 ±0.6) kpc 


ZlSM 


0.016 


5 


(0.012 ± 0.002) 



Table 5. Observed properties of the Milky Way compared with 
our predictions from a set of galaxies derived by applying the 
conditions presented in table [I] The references are as that table 
except: 4. |Binney fc Evans 20011 and 5. |Pagel 19981 The uncer- 
tainties in the fourth column are the ±1-(T limits of the distribu- 
tions. 

low and high circular velocity end, due to galaxies that have 
exhausted their gas supply. 

We leave a detailed consideration of galaxy sizes (bot- 
tom right panel) to the next section. 

8.4-1 Predictions for the Milky Way 

To provide a consistency check for our models, we will com- 
pare MW-like galaxies in our simulations with the known 
properties of our galaxy. On each panel of Fig. |5] we show 
a grey hexagon for the approximate location of the Milky 
Way itself. The values adopted for the Milky Way are pre- 
sented in tables 0] and |5] We also select MW candidates 
from our galaxies by applying the selection criteria shown 
in table El an d requiring that the galaxy have spiral (or Sb) 
morphology (defined in section F8.6H . Selection of this kind 
produces around six hundred MW galaxies in the simula- 
tion. The properties of these galaxies are shown in table 
where we quote the median and one standard deviation of 
the distributions, and compare to data for the Milky Way 
itself. We find a generally good agreement, in that applying 
the selection of table 01 produce galaxies that do look a lot 
like the Milky Way. For most of the properties considered, 
the Milky Way would fall within 1- or 2-a from the median. 

8.5 Galaxy sizes 

In the final panel of Fig. [H] we show the galaxy half-mass 
radius. This is obtained from a mass- weighted average of the 
disc and bulge half mass radii. The straight line is the radius 
of a galactic disc with the same velocity as its parent halo, if 
the radius is given by 1.68 Ai?20o/\/2, ie. equation 14.11 with 
A = 0.04. 

In Fig. [5] we examine the sizes of field disc galaxies (ie. 
galaxies which are alone in their halos) emerging from our 
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Figure 8. Key galaxy properties, as a function of characteristic velocity, which is a mass-weighted average of the disc circular velocity 
and the bulge dispersion velocity. The top-left panel shows the total baryonic mass of the galaxy, ie. stars plus cold gas. Top right is the 
metallicity of the interstellar medium, the ratio of the mass in heavy elements to the total gas mass. In the lower-left panel we show the 
star formation rate, and in the lower-right the galaxy radius, defined as the mass- weighted average of the disc, bulge and burst radii. For 
clarity, on this last figure we also plot the half mass radius of a disc forming in a halo with velocity V^ a i if its spin parameter were equal 
to 0.04 (equation 14.11 . In each panel, the grey hexagon shows the approximate location of the Milky Way, from data given in tables l4l 
and |5] 



model. We compare our predictions with the observed bivari- 
ate space density distribution of de Jong and Lacey ( 2000 ) 
derived from a sample containing more than a thousand 
field spiral galaxies. As stressed by these authors, the main 
advantage of using this data is that it should be quite in- 
sensitive to surface brightness incompleteness effects. Note 
that in the model we included all the spirals with bulge- 
to-total-light-ratio less than 0.33, i.e. morphological types 
later than Sab, whereas the data does not contain types 
later than Sd. Indeed these very late type galaxies cannot 
be excluded from the model sample as the only morpho- 



logical information we are able to compute is a bulge to 
disc ratio. We would therefore expect an overprediction of 
the total number of galaxies, especially in the faint magni- 
tude bins and for small disc sizes. However, we do not see 
this excess of small faint galaxies but rather an excess of 
bright small galaxies with respect to the data. Moreover, as 
pointed out by |de Jong fc Lacey (2000 | the model distribu- 
tion is (only slightly in our case) broader than the observed 
one. As we assume specific angular momentum conservation 
during disc formation, the obvious culprit is the too broad 
initial spin distribution of the dark matter halos. However, 
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Figure 9. Bivariate density distribution of field spiral galaxies 
as a function of effective disc radii in different I-band absolute 
magnitude bins. Magnitude cuts are indicated in the top corner 
of each panel. The solid histogram represents results from the 
hybrid model, and the solid line the best fit to the data of de 
Jong and Lacey (20001 (their cq. 7). The horizontal dashed lines 
in each panel represent our volume resolution. 



we would like to draw the attention of the reader to the 
fact that unlike the CLBF model which overpredicts quite a 
lot the density of small discs, we seem to predict about the 
right number of small-sized bright galaxies. Although our 
mass resolution biases us against small-sized galaxies, we 
estimate (see table [^J that our sample is almost complete 
for galaxies with absolute magnitudes brighter than -20.9 in 
the I band. Therefore the results displayed in the first three 
panels of figure [5] should not be too sensitive to an increase 
in resolution. 

Once again we emphasize that our recipe for disc 
sizes depends crucially on the assumption of specific an- 
gular momentum conservation, and it has been shown 
(eg. Navarro, Frcnk, & White 1995a) that disspative ef- 
fects can play a part in destroying some of the an- 
gular momentum of the gas, producing smaller discs. 



Dominguez-Tenreiro, Tissera, & Saiz (1998 i show that this 
loss is less important for gas cooling onto 'evolved' galaxies 
with stellar bulges, so these effects are likely to preferentially 
reduce the size of the smaller, less resolved galaxies in our 
model. 



8.6 Morphologies 

We predict galaxy morphologies based on the ratio of 
B-band luminosities of the disc and the bulge compo- 
nents. It has been found ISimien fc de Vau couleurs 198611 
that this ratio correlates well with Hubble type. We de- 
fine a morphology index, / = exp(— Lb/Lq), such that 
a pure bulge has value 0, a pure disc 1. Following 
|Battgh, Cole, fc Frenk (1996] , we translate this index into 



a morphology by assuming ellipticals have I < 0.219, SOs 
have 0.219 < I < 0.507, and spirals have J > 0.507. 

We first consider the morphological mix of galaxies in 
our simulation. Using the APM Bright Galaxy Catalogue 
(BGC, |Loveday 19961 , |Baugh, Cole, fc Frenk (T996| > point 
out that in the local Universe, the ratio E : SO : SP + Irr = 
13 : 20 : 67. In fact, one cannot directly compare this to 
the ratio arising in our simulations, since the catalogue is a 
flux-limited sample of the Universe, whereas our samples are 
effectively volume limited, containing an entire simulation 
cube at z = 0. If there is any correlation between B-band 
luminosity and morphological type, some types may be seen 
over a larger volume in the BGC than others. 

For a better comparison, we take the Stromlo-APM 
redshift survey fLov eday et al. 1996^ and construct from it 
a volume limited sample of galaxies containing all those 
with well-measured morphology, down to the absolute mag- 
nitude limit of our models, Mb < —18.8. We obtain 
E : SO : SP + Irr = 13 : 11 : 76 (the sample contains just 120 
galaxies, so the Poisson errors are ±3 : 3 : 8). Thus it ap- 
pears that, since SO's are intrinsically brighter objects than 
spirals, the ratio from the BGC overestimates their abun- 
dance when this selection is made. We compare to the ratio 
found for our models with a similar magnitude limit. We 
find 17 : 16 : 67 for a sample of twenty thousand galaxies, 
with a total number density around twenty per cent higher 
than that measured in the volume limited Stromlo-APM 
sample. Thus it seems that for these galaxies we slightly 
over-represent Es and SOs and under-represent spirals, al- 
though it must be remembered that the statistics from the 
real data are poor. 

If we apply a more conservative limit, going 0.75 mag- 
nitudes brighter (ie. a factor two in luminosity), we find 
the Stromlo-APM data give 14 : 13 : 73 (250 galaxies in 
the sample, since the volume probed has increased, which 
implies Poisson errors 2:2:5), and 21 : 17 : 61 from our 
models. Again our number density is around twenty per cent 
higher than in the real data. In other words, as we increase 
the absolute magnitude of the selected objects, the quanti- 
tative comparison with the data remains quite stable. This 
effect is unsurprising, since in section |8~T1 we stated that our 
magnitude for completion in the B-band is around -18.9. 

It is worth noting that, with more complete surveys such 
as 2dF and SDSS, volume limited samples of the local Uni- 
verse with much better statistics will be available, and these 
numbers can be checked more accurately. As it is, we are 
thus broadly confident that our methods produce approxi- 
mately the correct morphological mix at low-redshift. How- 
ever, it must be appreciated that we are comparing a simple 
recipe for determining morphologies with what is observa- 
tionally a sophisticated and subjective process (ie. assuming 
a one-to-one mapping between bulge-to-disc ratios and Hub- 
ble type), and one should not place too much emphasis on 
these results. 



8.7 Colours 

We next consider the dependence of colour on galaxy mor- 
phology. Fig. HOI presents B — V colours as a function of mor- 
phology for our galaxies. The solid line represents data taken 
from lButa et al. 19941 Our spiral galaxies closely follow the 
data both in average value and scatter, but more important 



22 Hatton et al. 




exp(-B/D) B 

Figure 10. B — V colours as a function of morphology for our 
galaxies brighter (than -18.8 in the B-band) seen face-on. The 
solid triangles show the median of the distribution, with the 
dashed error bars giving the 1-cr range. The solid line shows the 
measurements of Buta et al. 1994 (converted to a bulge-to-disc 
ratio index according to Simien &; dc Vaucoulcurs 1986 I with the 
observed l-cr dispersion (dot-dashed lines). The short-dash verti- 
cal lines show where the divisions are drawn between elliptical, 
lenticular and spiral morphologies. 

is the dependence of colour on morphological type for early 
type galaxies which is not as strong as observed, leading to 
elliptical/SO galaxies which are too blue with too strongly 
scattered colors, independently of the exact criterion used to 
define these morphological types. We will come back to this 
important issue in paper II, but we mention that this is due 
to recent overcooling leading to late star formation in early 
type galaxies. Indeed, if we remove ellipticals and SOs with 
more than 1 % of gas in mass still present at z — our colors 
(and dispersions) are in good agreement with the observa- 
tions, although more than half of our early-type galaxies 
drop out of the sample. Finally, we note that even though 
this discrepancy is problematic, our median colors are still 
within 1-2 a of the observations, whatever the morpholog- 
ical type considered. Furthermore, RAM pressure stripping 
could help redden the colours of SOs and cluster ellipticals. 

8.8 Luminosity functions 

The luminosity function, ie. the number density of objects 
as a function of magnitude, is a commonly measured obser- 
vational statistic. In Fig. llll we show the results for the lumi- 
nosity functions arising from our models in different wave- 
bands. We compare these with data coming from several dif- 
ferent sources. For the optical and near infrared (.B-band and 
K-band) we use recent data from the 2dF galaxy redshift 
survey llCross et al. 20011 for the B-band, lUole et al. 200T1 
for the K-band). For the IRAS wavebands (25 and 100 /kid.), 
we compare with results from |Soifer fc Neugebauer (199l[ . 
As noted in section 18.11 the finite mass resolution of 



our simulations produces a lower bound to the galaxy lu- 
minosity to which we are complete in any given waveband. 
Although we have galaxies fainter than this limit, there will 
be galaxies of similar luminosity in halos that we do not 
resolve, and hence we expect to see some depletion of the 
counts lower than this limit. We plot the magnitude limit 
corresponding to the mass limit with a vertical dashed line 
in each of the four plots in Fig. 1111 The resolution effect, ie. 
the depletion of counts faintwards of this line, is apparent in 
all four wavebands. For the optical and near infrared, this 
depletion seems to occur quite close to the resolution limit, 
starting 0.5 - 1.0 mag brighter. This suggests that this is 
an accurate estimate of the resolution limit in these bands. 
In the IRAS bands, the turnover in the counts occurs at a 
luminosity of around 0.5 dex (roughly 1 mag) brigther than 
the formal limit. 

On the whole, the plots show quite good agreement 
with observed galaxy number densities over two orders-of- 
magnitude in wavelength. In the B-band, the counts agree 
reasonably well with the shape and normalization of the 2dF 
data. This is a significant improvement on the hybrid model 
results of KCDW who found a fairly fiat, power-law slope for 
the B-band luminosity function in a ACDM Universe. Our 
results are similar to those of CLBF, although they are able 
to get a slightly better amplitude of the luminosity function 
because they have a free parameter to control the fraction 
of stars forming brown dwarfs (which of course do not con- 
tribute to the luminosity function), that they fix using the 
B-band number density at Mbj — 51og 10 h = —19.8. 

In the AT-band, we get a good (although on the high 
side) overall number density, and fit the data quite well at 
the 'knee' of the luminosity function, representing most of 
the old stellar mass in the Universe. Moreover, we note that 
our overprediction of intermediate luminosity galaxies with 
respect to the 2dF data might not be so much of a problem 
since a more recent survey with similar number statistics 
( Hua ng et af. 2003| l seems to also find a higher number den- 
sity of these galaxies. Our luminosity function has an overall 
appearance which is similar to the Schechter function that 
well describes the data ie. we reproduce the exponential cut- 
off at the bright end that is present in the data. This be- 
haviour is due both to the prescription we employ to solve 
the overcooling problem (see section 13.31 and the more ef- 
ficient feedback we implemented (see section T4. 21 . KCDW, 
with their different prescriptions for both these effects, ob- 
tain a power-law shaped K-band luminosity function which 
overestimates the number of bright objects. CLBF obtain a 
better fit to the data, which is insensitive to dust emission 
and absorption but once again relies on the free parameter to 
control the number of brown dwarfs formed. The conclusion 
is that overpredicting the bright end of the if-band luminos- 
ity function seems a rather natural feature of semi-analytic 
models and one has to invoke at least a rather strong feed- 
back (from supernovae and/or AGNs) to get rid of stars in 
the the brightest galaxies. 

For the IRAS number densities, we again find fairly 
good agreement with available data except at the faint end 
where resolution effects dominate. In both wavebands the 
agreement is of excellent quality over the whole range of 
luminosities above the formal resolution limit. 

On the whole, we find the broad agreement with obser- 
vations over several orders-of-magnitude in galaxy luminos- 
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Figure 11. Galaxy luminosity functions. Top left: the B-band luminosity function (points with error bars), with data from the 2dF 
galaxy redshift survey (Cross ct al. 20011 (connected, lighter curve with error bars). Top right: the i\-band luminosity function, with 
data from 2dF ICole et al. 200U . The error bars are Poisson errors, and the magnitudes are given in the AB system. Bottom left: IRAS 
25 (im luminosity function. Bottom right: IRAS 100 fim luminosity function. fflylSdata comes from |So'i fcr & Ncugebau er 199l] In each 
panel, the dashed vertical line is the formal resolution limit defined in section f8.il and the dashed horizontal line is the volume limit, (ie. 
one galaxy per magnitude bin). 



ity and two orders-of-magnitude in wavelength, a good con- 
firmation of our models both for the semi-analytic recipes of 
the physics of galaxy formation, the determination of spec- 
tra, and the modelling of extinction effects. 



8.9 Optical— IR colour magnitude relation 

In Fig. 1121 we examine the behaviour of the optical-IR 
'colour', ie. the difference in luminosity between these wave- 
length ranges, as a function of galaxy luminosity. We include 
all galaxies with baryonic masses greater than the formal 
mass limit, as defined in section IST1 Soifer et al. (1987|l find 



that the colours of bright galaxies in the IRAS sample are 
closely related to their total IR luminosities, but have no 
correlation with the optical luminosity, implying that galax- 
ies have similar optical properties no matter what their IR 
emission is. It can be seen from Fig. EU that our galaxies 
follow a similar trend. On the left hand side we show the 
dependence of the ratio of bolometric IR luminosity to B- 
band luminosity against the IR emission. For low to moder- 
ate luminosities, this relationship is close to linear. A linear 
regression reveals that the slope of this correlation is 0.70. 
Fitting a linear relationship "by eye" to the data of Soifer 
et al. (their fig. 5a) yields a slope of 0.75, so the results 
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Figure 12. The relationship between the optical-IR colour, defined as the ratio of infrared to blue luminosity, and the absolute bolomctric 
luminosity in the infrared (left) and S-band (right). All galaxies with baryonic masses greater than the formal mass limit are plotted 
(section 18.11 . Solid lines represent lines excluding about 1 a of the data plotted in figs 5a and 5b of Soifcr ct al. 1987 



are quite comparable. Our distribution is somewhat skewed, 
with a population of galaxies of quite high IR flux densi- 
ties and extremely red colours, departing from the slope 
that fits the rest of the galaxies quite well. Since galaxies 
in Soifer's plot need both IR and B-band luminosities to 
be measurable, galaxies that are intrinsically faint in the 
B-band are discriminated against in this sample. It is thus 
likely that many of the galaxies in the upper part of our dis- 
tribution would not be included when this sort of selection 
is applied. This is especially clear in the right hand panel 
of Fig. 1121 where we plot the colour against the B-band lu- 
minosity. No correlation is visible, except for the handful of 
extremely red galaxies, which all have B-band luminosities 
less than 3 x 10 9 L Q . Note that in lSoifer et al. 19871 (their 
fig. 5b) there are very few galaxies indeed with B-band lu- 
minosities this low. Furthermore, the IRAS sample is flux 
limited whereas ours is volume limited, and therefore the 
observations should be biased towards IR bright and dis- 
tant objects, which partly explains why our population of 
galaxies seem to populate the low-IR part of both diagrams. 
The other, and more fundamental reason is that although 
we only plot galaxies above the formal mass resolution limit, 
we do not satisfactorily resolve the star formation history of 
galaxies just above this threshold, which obviously leads to 
bluer colors and lower IR emission because according to the 
simulation, most of them have formed very recently. 

Bearing these limitations in mind, we conclude that our 
model is fairly successful at reproducing the observed cor- 
relations between luminosity and colour in the infrared and 
blue wavebands. 

8.10 Tully-Fisher relation 

There is a well-established empirical relationship be- 
tween the luminosity and rotation velocity of galactic 
discs. This effect is known as the Tully-Fisher relation 
HTully fc Fisher 1977) . In the top left hand corner of Fig. [151 



we present a scatter plot for the /-band magnitude against 
circular velocity for all spiral galaxies in our simulation. The 
dashed line on the plot represents the best fit to the I- 
band Tully-Fisher relation found by |Giovanelli et al. (1997| l , 
where we have assumed that the 21cm line width is simply 
twice the galaxy circular velocity. We plot the relation only 
for the range of velocities over which Giovanelli et al. have 
data. We note that the value of the Hubble constant derived 
from this dataset is h — 0.69 ±0.05, which is consistent with 
our value of 2/3. It is clear from this plot that the major- 
ity of galaxies in our sample do indeed fall in a single locus 
on these axes, although there is significant scatter in the 
relation and a population of galaxies with rather lower lu- 
minosities for a given circular velocity. The main locus itself 
is close to the relation of |Giovanelli et al. (1997^ although 
our galaxies appear to be somewhat fainter, and the slope 
of the correlation is shallower. 

We proceed to perform a more detailed analysis by split- 
ting our galaxy sample into three subsamples. These are: 

• 'field' galaxies, that are the sole occupants of their dark- 
matter halos, 

• 'cD' galaxies, which are defined as central galaxies in 
halos with more than one occupant, 

• 'satellites', all non-central inhabitants of multi-member 
halos. 

In the other three panels of Fig. 1131 we show the Tully- 
Fisher relation for each of these subsamples. In each case, 
we plot the line connecting the median of the binned distri- 
bution, and add error bars to show the ±1-<t range. We again 
compare to the line from |Giovanelli et al. (1997| l. It can be 
seen from the plots that our cD galaxies follow the observed 
Tully-Fisher relation quite well, although these galaxies are 
slightly too faint (or rotate too fast). 

It is clear that the observational Tully-Fisher relation is 
not followed nearly so well by our field and satellite galaxies 
(the solid lines in the panels of Fig. I13H . The median mag- 
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Figure 13. The Tully-Fisher relation in the /-band for all our spiral galaxies above the resolution limit (top left). The other three panels 
show the same data split into subsamples by galaxy environment, with the median and rtl-<7 scatter. In each panel, the dashed diagonal 
line is the observational Tully-Fisher relation derived by |Giovanelli et al. (1997) . We also show the distribution both before (dashed line) 
and after (grey solid line) making the gas column density selection of equation 18. T1 



nitude as a function of circular velocity is well below the 
observed trend, and the 1-cr scatter is very large compared 
to that for the cD galaxies. Also the satellite galaxies are 
largely responsible for the large scatter in the plot in the top 
left. This is especially interesting because the sample used 
bv lOio vanclli et al. 1997 to measure this relation is actually 
a sample of spiral galaxies in 24 galaxy clusters, ie. corre- 
sponding to our satellites. There are two key reasons for this 
discrepancy: first is that satellites or galaxies with massive 
bulges are no longer being 'fed' with hot gas cooling onto 
the discs, with the result that they use up their store of cold 
gas, and consequently have low star formation rates and lu- 
minosities. In other words, galaxies that have been satellites 
for different lengths of time can therefore have quite differ- 



ent gas fractions, and this is responsible for the large scatter 
in the relation. Second, and more fundamental is that the 
circular velocity of a galaxy results from the complex inter- 
action between baryons and dark matter, and as a result our 
simple treatment of the disc dynamics is probably too crude, 
especially for the fast rotating (massive) systems which are 
also likely to have a complex merging history. 

As pointed out by CLBF, the data in fact are subject 
to a strong selection effect in that galaxy rotation velocities 
are measured using HI or optical emission lines, and a sig- 
nificant column density of interstellar gas must be present 
in order to make those measurements. In fact, many of our 
satellite spirals are gas-poor, and so would not be included 
in the IGiovanelli et al.l dataset. To briefly demonstrate this 
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selection effect, we applying a cut-off mimicking that which 
must exist in the data, 

M cold /Tvr 2 1/2 > f0 7 M Q kpc- 2 , (8.1) 

which corresponds to a column density of approximately 
iV H = 1.25 10 21 atoms cm . This selects only the satellites 
that still have a significant amount of gas, and the effect 
is clearly visible on the Tully-Fisher relation (solid line in 
Fig.[inj. The median now has similar behaviour to that of 
the field galaxies, although the satellites are still between 
0.25 and 0.5 magnitude fainter than field galaxies depending 
on circular velocity. The scatter has been vastly diminished, 
although it is still larger than that of the field galaxies. This 
also reduces the number of satellite spirals by about a fac- 
tor of 2.5, in keeping with a number of observational results 
(eg. ISolanes et al. 2001ft demonstrating that spirals in clus- 
ter environments are often gas-depleted. 

We can also compare the scatter in the Tully-Fisher 
relation with the intrinsic scatter estimated from the ob- 
servations of IGiovanelli et aTl We compare the scatter for 
a MW-type galaxy with circular velocity 220km s~\ using 
their equation 11. Their prediction is a scatter of £i nt = 0.22 
magnitudes. We find a slightly higher scatter for our field 
galaxies, by interpolating a spline fit to the error bars shown 
in Fig. 1131 eint = 0.33. The 'cD' galaxies have a similar scat- 
ter, tint = 0.24, and the cluster spirals have much larger 
scatter than the observations predict, etot = 0.97 for the 
whole sample and ejnt = 0.49 after applying the column 
density selection of equation 18.11 Whilst it is encouraging 
that our cD galaxies have similar scatter than the data, it 
seems that even with selection, our model predicts a wider 
variety of spirals in cluster environments than are observed. 
These results imply that, in reality, either gas cooling plays 
the dominant role in determining the properties of spirals in 
both cluster and field environments, and our suppression of 
cooling onto satellite and massive field galaxies is responsible 
for their variety, or our estimates of galaxy velocities are sys- 
tematically offset towards higher values and more dispersed 
than in reality. 

8.11 Faber— Jackson relation 

Analogous to the Tully-Fisher relation for spiral galaxies, 
observations show a correlation between the luminosities 
and dispersion velocities for elliptical and SO galaxies, known 
as the Faber-Jackson relation IIFaber & ; Jackso n" 1976t . In 
Fig. ll4l we show the Faber-Jackson relation for our galaxies, 
and compare with that found by |Forbes fc Ponman ( 1999| > 
for cluster ellipticals (shown by the solid line). We mimic 
the observational selection by plotting only elliptical and SO 
cluster (i.e. halos containing more than 10 galaxies) galaxies 
but note that the match is qualitatively similar if we use all 
the galaxies we classify as elliptical or SOs. 

It will be seen that the galaxies sit around a line which 
is slightly fainter (or at a slightly higher velocity dispersion) 
than that of Forbes & Ponman, once again consistent with 
too high velocity dispersion estimates. However we consider 
this a fairly good imitation of the observed Faber-Jackson 
relation for our galaxies because of the scatter in the rela- 
tion, which is shown as an inset on this figure. This has been 
obtained from the differences between the straight line and 
the individual velocities as a function of B-band magnitude. 



2.6 




1.8 



-18 -19 20 -21 -22 -23 



Figure 14. The Faber-Jackson relation for our elliptical and 
SO cluster galaxies. We only plot galaxies above our resolu- 
tion limit Ms = -18.8. The solid line shows the fit found by 
|Forbes &: Ponman (1999} for local ellipticals. The inset shows the 
scatter of our galaxies from the observed relation. 

We have shaded the central area containing 68.3 per cent of 
the data, and we can thus read off the 1-a limit as +0.12 and 
—0.15 dex on each side of the observed relation, correspond- 
ing to approximately ±0.13 dex of scatter if we were to define 
a Faber-Jackson relation directly from our model galaxies. 
Estimates of the scatter by Forbes & Ponman yields a 1-a 
limit of ±0.10 dex, so we are confident that even though our 
relation looks slightly offset from the data our models still 
are in fairly good agreement with observations. 

8.12 Fundamental plane 

It has been found observationally that the scatter in the 
Faber-Jackson relation is rather larger than what would be 
expected simply from measurement errors, leading to the 
idea | |Djorgovski fc Davis 1987| IDressler et al. 19871 that 
there is a further variable involved. The so-called Funda- 
mental Plane relates the galaxy luminosity, velocity disper- 
sion, and radius, and is found to apply to elliptical galaxies, 
and, to some extent, the bulges of spirals as well. 

Following a similar procedure to that employed by ob- 
servers, we assume that the plane is of the form: 

log R e = a log a ± bfi c ± c (8.2) 

In this equation, R e is the effective radius of the spheroid, 
defined as the radius that contains one-half the galaxies 
light when viewed in projection (for our bulges this is _R C = 
1.815 ra, see Hcrnquist 1990). a is the one-dimensional ve- 
locity dispersion, and p e is the average surface brightness 
within R c . 

We determine the coefficients (a, b, c) in eauation l8.2l bv 
minimizing the total absolute deviation of all our elliptical 
and SO galaxies above the formal resolution limit from this 
plane. In order to do this, the variables are first transformed 
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Figure 15. The fundamental plane for our elliptical and lenticular galaxies, in the B-band (left) and /-band (right), above our resolution 
limit. The coefficients on the y-axis are found by minimizing the absolute deviation from the line y = x, as explained in the text. The 
insets show the distribution of the residuals about the plane. 



into a 'normalized' co-ordinate system by subtracting off 
their mean value and dividing by the standard deviation, ie. 



0i = (<r — a)/ a 



(8-3) 



We show our Fundamental Plane relation in Fig. 1151 in 
both the B-band and the /-band. There is clear evidence for 
a fairly tight correlation of these three bulge properties, in 
both wave-bands, with the relation given by: 



(a, b) = 



(1.47,0.38) B-band 
(1.48,0.36) /-band 



(8.4) 



For a sample of 74 early-type galaxies towards the Coma 
cluster, |Scodeggio et al. (1998} find that the fundamental 
plane that minimizes the absolute deviation is given by 
(their table 6): 



(a,b) 



(1.40, 0.35) B-band 
(1.70,0.33) /-band 



(8.5) 



These results for the coefficients (a, b) are in good agree- 
ment with those we find (and reasonable agreement with 
other observations, eg. |Pahre fc Djorgovski T997J. It will be 
noted that we also succeed in replicating a trend that exist 
in the data as well, namely that the parameter b is quite 
independent of the waveband. The steepening of the slope 
a of the dependence on the velocity dispersion as a function 
of increasing wavelength, is not so marked in our simulated 
galaxies. 

In the insets to the two panels in Fig. 1151 we show the 
distributions of the residuals about the plane, ie. the differ- 
ences between the actual log R c values and the prediction 
from the derived FP relation. The width of these distribu- 
tions are approximately ±0.15 and ±0.12 for the B-band 
and /-band respectively. An estimate from Scodeggio et al. 
fig. 7 shows the 1-a scatter in the data to be approximately 
±0.10 dex in both the B-band and the /-band. Thus we not 
only produce fundamental planes with a similar tilt to that 
observed, but also obtain a good match for the residuals. 



This also suggests that much of the scatter in the observed 
FP relation is intrinsic, due to the different merging histo- 
ries of the objects examined, rather than being caused by 
measurement errors. 

This is the first time that semi-analytic models have 
been pushed as far as looking at the fundamental plane, 
and we find the results an encouraging vindication of our 
methodology. 



9 RESOLUTION TESTS 

The issue of resolution enters our models in several ways, 
and it is important to check that we understand how reso- 
lution affects the results presented here. In this section, we 
will endeavour to show that our results are quite stable to an 
increase in resolution. This will generally be done by com- 
paring with a simulation of slightly worse resolution, and 
showing that the results have already converged. 

In the following sections, we consider the effects of fi- 
nite force, time, and mass resolution. We will use the ab- 
breviations HTR and LTR to stand for high and low time 
resolution, and HMR and LMR for high and low mass res- 
olution. Thus the HTR (or HMR) is generally the original, 
high resolution simulation. 



9.1 Spatial resolution 

Initially, in the iV-body simulations, there is a spatial res- 
olution set by the softening length used in computing the 
smoothed gravitational potential. In our code, this softening 
length is set at one twentieth the mean interparticle sepa- 
ration. For our ACDM simulation, at redshift zero, this is 
around 30kpc. We identify halos with a density contrast of 
200 relative of the critical density, which thus have an aver- 
age particle-particle separation length of the mean separa- 
tion divided by (200/fit) 1//3 , where fi t is the ratio of matter 
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density to critical density as a function of time. The mean 
distance between particles inside halos at z — is thus twice 
the softening length. Improved resolution may then have an 
effect on the central, dense regions of halos, where the sep- 
aration length is smaller than the softening length, but is 
unlikely to influence the virial mass itself. Thus, attempts 
to fit density profiles to these inner regions would be af- 
fected by the softening, but we do not do this, measuring 
only virial masses and energies, and assuming an isothermal 
profile for all our halos. The spatial resolution limit set by 
gravitational softening is thus not a relevant concern for this 
work. 



9.2 Time resolution 

To calculate the particle trajectories in the A^-body simula- 
tions, we use around twenty thousand timesteps (table 1), 
equally spaced in expansion factor from z ~ 35 to the 
present day. To construct the halo merging tree, a subsam- 
ple of these steps is used, spaced logarithmically in expan- 
sion factor. This produces seventy timesteps between the 
present day and the first identification of a virialized object 
with twenty particles, at redshift z w 10. The spacing be- 
tween the penultimate and the final timesteps corresponds 
to sixty of the simulation timesteps. It is thus the choice 
of output times for constructing the merging tree that re- 
ally dictates time resolution. This especially has an effect 
on the baryon cooling tree described in section YTM since we 
will only be interested in binary mergers in this section. If 
a massive halo has accreted two smaller halos in the course 
of one timestep, we only follow the galaxies in the larger of 
the two halos. The occurrence of more-than-binary mergers 
obviously increases with the length of the timestep, so more 
galaxies are lost with worse resolution. 

To check our sensitivity to this timestep resolution, we 
compare results with those from a tree made using only one 
half the number of simulation outputs (ie. with a factor two 
worse time resolution). There is a one-to-one correspondence 
between halos (but not galaxies) in the final timestep of the 
two trees, so we can look to see how resolution has had an 
effect on the properties of individual halos. 

Fig. Ilbl shows a halo-by-halo comparison for several key 
halo properties. The solid line shows the number of objects 
with a given fractional difference between their cold gas 
masses in the full simulation and in the low time resolution 
version. It can be seen that there is a very small systematic 
offset between halos in the two simulations, with halos in 
the HTR having slightly more cold gas, and there is some 
scatter (the curve has a FWHM of the order ±20 per cent). 

For hot gas, the picture is rather different. Although 
there is a sharp peak at zero difference, this is due to halos 
that have just been discovered in the final timestep, and 
so have not cooled at all in either version. For the rest of 
the halos, the peak in the distribution occurs when halos in 
the HTR run have around forty per cent more hot gas than 
those in the LTR. However, the total amounts of hot gas in 
the two simulations are similar, with the HTR run having 
around eight per cent more gas than the LTR. This implies 
that the halos containing large amounts of the hot gas are 
the same in the two simulations, but that a large number of 
halos containing small amounts of gas are discrepant, and 
their precise hot gas mass is sensitive to the time resolution 
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Figure 16. The difference between several halo properties (com- 
puted on an individual basis) between halos in the full resolution 
(HTR) run and halos in the run only using half the number of 
timesteps (LTR), as explained in section flJ. 21 
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Figure 17. Differences in B-band magnitudes for individual halos 
in the full and half time resolution runs. The upper (solid) line is 
the magnitude difference of the brightest halo member (BCG) for 
halos with mass less than 3 10 12 Mq. The middle (dotted) line is 
the same statistic for the other, more massive halos. The lower 
(dashed) line is the difference in satellite magnitude (ie. summed 
over all galaxies in the halo excluding the brightest one) for each 
halo having more than one member galaxy. A vertical line has 
been plotted at zero difference to guide the eye. 
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chosen. It is interesting to note that this effect does not have 
a noticeable impact on the metallicity of the hot gas, shown 
by the short-dashed line in Fig. 1161 the metallicity difference 
is peaked close to zero, showing that whatever mechanism is 
responsible for losing gas in the LTR run also loses a similar 
fraction of metals. With the long-dashed line, we show the 
difference in total stellar mass for each halo. Happily, this 
is peaked close to zero, so it seems that time resolution has 
little effect on the total stellar masses. 

In Fig. 1171 we present three curves comparing luminosi- 
ties between the two runs. The upper curve is the B-band 
magnitude difference between the brightest galaxy in the 
halo (BCG) in the two different runs, for halos of virial mass 
less than 3 1O 12 M0, although the majority of halos in this 
mass range in fact only contain a single galaxy. It can be 
seen that this is peaked to the left of zero, implying brighter 
BCG's (by around 0.2 magnitudes) in the HTR run. This is 
despite the fact that we know from above there is no offset 
in the total mass of stars between the two runs. Since it is 
the B-band magnitudes we are comparing, which are domi- 
nated by contributions from the youngest stars, this implies 
that the long timestep in the LTR smoothes out the star 
formation over too long a period. The second (dotted) curve 
compares the BCG's in the other, more massive halos. The 
offset is reduced (~ 0.1 magnitudes) so this effect is less im- 
portant for massive halos. In the lower curve, we show the 
difference in total halo magnitude when summed over all 
galaxies in the halo except the brightest one. Again, there is 
more luminosity in the HTR run. In fact, the final timestep 
in the HTR run contains a total of 6355 satellite galaxies, 
as opposed to 5096 in the LTR. This is because galaxies are 
lost due to more-than-binary merging in the LTR run, as 
explained in at the beginning of this section. This twenty 
per cent difference in galaxy numbers is largely responsible 
for the ten per cent difference in satellite luminosity in this 
plot. 

In summary, there would have been little difference had 
we run our simulation with worse time resolution, except for 
the depletion of low-mass satellites and slightly fainter field 
galaxies. Since these effects are all of order 0.1-0.2 mag- 
nitudes, which is quite small compared to the systematic 
errors due to uncertainties in the modelling techniques, we 
conclude that our time resolution is close to convergence. 
Full convergence could be achieved with better time reso- 
lution, though we have chosen in the previous parts of this 
paper to gain improvement by allowing the code to cope 
with multiple merging of halos. 

9.3 Mass resolution 

There are two ways in which mass resolution enters our re- 
sults. Initially, in the A-body simulations, there is a res- 
olution limit set by the finite mass of each particle. Sub- 
sequently, when identifying dark matter halos, there is a 
minimum number of particles needed. We will consider the 
impact of both these effects. 

9.3.1 N-body resolution 

To investigate the effect of A-body resolution, one would ide- 
ally run two simulations of different resolution and compare 



the results. Although we have runs of 64 3 and 128 3 particles 
for the same box size and initial conditions as the standard 
(256 3 ) simulation, there is very little overlap in the halo 
mass function between them, making a direct comparison 
rather unproductive. Instead of comparing the simulations, 
we degrade the resolution of the HMR run by selecting, at 
random, one half the particles, assigning them double their 
original mass. We then run the halo-finding, tree-building 
and galaxy formation codes on this LMR (degraded) sim- 
ulation. In Fig. 1181 we show how the luminosity functions 
(in the B-band and IRAS 60 /im filter) of these simulations 
differ. For faint galaxies, the LMR run underestimates the 
number density by a factor ~ 2.5. For bright galaxies, the 
estimates converge. It is the behaviour between these two 
regimes that is especially interesting. The vertical line in the 
plot is the formal resolution limit {Mb — 51og 10 h < —18.8 
or M6o M m — 51og 10 /i < 22.9, see section of the LMR 
run, and it will be seen that it is right in the middle of the 
transition between these two regimes, and that the LMR 
simulation does not fully converge with the HMR until ap- 
proximately three-quarters of a magnitude brighter than this 
limit. As stated in section we expect to have problems 
for slightly brighter galaxies than this limit, since we can- 
not accurately follow the merging histories of the lightest 
objects. 

Assuming that this behaviour is similar whatever the 
resolution limit, we have therefore shown that our galaxy 
'selection' is complete for galaxies around three-quarters of a 
magnitude brighter than the formal resolution limit. Hence, 
our results are reliable for galaxies down to one magnitude 
fainter than L*. 

This is not a direct test of the effects of A-body res- 
olution, since we have compared with a degraded version 
of the same simulation rather than a simulation actually 
run with fewer particles. In order to check the validity of 
this degrading approach, we sample the 256 3 simulation at 
a much lower rate, selecting one particle in 4 3 , thus giving it 
the same resolution as the 64 3 simulation. Making a similar 
comparison, we find no significant differences in the lumi- 
nosity functions. Having verified this approach with a factor 
of 64, it is reasonable to expect that sampling with a factor 
two will produce a good imitation of a simulation run with 
half the mass resolution. 

9.3.2 Minimum particle number 

Increasing the minimum particle number for halos will mean 
that the building blocks of our galaxy formation process 
are more massive, so the effect should be similar to that of 
degrading the mass resolution of the simulation. However, 
the halos will be better resolved in this new case, so there 
may be some differences. To investigate this, we compare the 
LMR (degraded) simulation of sect ion 19 . 3 . 1 1 wit h the HMR 
simulation when the minimum particle number is increased 
from twenty to forty. In Fig. 1191 we compare the B-band 
luminosity functions of galaxies in these two runs. As can 
be seen, there is no significant difference between the two. 

The overall conclusion is that the limit of our resolution, 
in terms of noticeable effects in the local Universe, depends 
only on the mass of the smallest halos identified, and is not 
individually sensitive to the actual numerical resolution of 
the simulation or the number of particles per halo. 
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Figure 18. Differences in the B-band (left panel) and IRAS 60 (im (right panel) luminosity functions for the full resolution and half 
resolution runs. The vertical line is plotted at the estimated resolution limit of the LMR simulation (approximately 0.75 magnitudes 
brighter than values for the HMR simulation quoted in table IS! . 




Figure 19. As Fig. 1181 but this time we compare the luminosity 
function from a run with a minimum particle number per halo of 
forty, with that of the half-resolution sampled run. 



This also suggests that it would be valid to improve the 
resolution of these runs by extending the treatment to halos 
with less than twenty particles in them, for example going 
down to as few as seven, the limit quoted by KCDW for 
dynamical stability. This alone should result in roughly a 
one magnitude improvement in resolution. 



10 DISCUSSION 
10.1 Summary 

In this paper, we have presented our A^body/semi-analytic 
hybrid model for studying galaxy formation, and demon- 
strated its application to a reference model within a ACDM 
cosmology. Here, we summarize the assumptions and weak- 
nesses of the model, emphasize some of its results, and com- 
ment on several aspects of the programme. 



10.1.1 Implementing new recipes 

At root, our approach is entirely independent of pre-existing 
work as we use different N-hoAj codes and different methods 
for tracing the history of dark matter halos. However, most 
of the recipes used to model physical processes (eg. gas cool- 
ing, star formation, feedback, dynamical friction) are similar 
to those that have appeared elsewhere. We include a number 
of innovations that have not previously been used in hybrid 
models: 

• We have measured the halo properties directly from 
groups identified in the TV-body simulation, rather than as- 
suming distributions eg. for the spin parameter. 

• We have noticed that the friend-of-friend group finder 
identifies a number of (mainly low-mass) halos that are not 
bound systems. We have dealt with these objects by sup- 
pressing cooling in them. If they are simply random collec- 
tions of particles, misidentified, they will thus not affect our 
results. 

• In general, we have endeavoured to use 'natural', con- 
tinuous models as opposed to the all-or-nothing methods 
applied by other workers. Examples include the model for 
halo reservoirs, where metals and hot gas ejected from the 
halo are retained in the local environment and gradually re- 
accreted, and the recipe for merging, where the effects are a 
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smooth function of the mass ratio of the progenitors, rather 
than the discontinuous step function. 

• We have modelled disc instability as another possible 
mechanism for the formation of bulges. 

• We allow for mergers between satellite galaxies in clus- 
ters, and gradual tidal stripping of dark matter cores of sub- 
halos. These follow schemes that have been applied to pure 
semi-analytic models by SP99. 

• We use the detailed chemical and stellar evolution mod- 
els described in DGS, giving us self-consistent, metallicity- 
dependent spectra and cooling. We do not make the common 
assumption of instantaneous recycling of metals. 

A few key weaknesses of our recipes present themselves: 

• The process used to subdue cooling in massive halos, 
and in small halos after recent mergers, is critical to avoid 
overcooling on to the central object, and hence overly bright 
central galaxies. Our approach has been to treat this prob- 
lem phenomenologically but to point out in details where it 
affects our final results. The main reason behind this choice 
being that whilst an approach like that of CLBF can be 
qualitatively justified, and is perhaps more reasonable than 
the simple approximation of forgetting about cooling in ha- 
los above a certain circular velocity (KCDW), there is no 
quantitative or observational justification for it. 

• The treatment of the interactions between dark mat- 
ter and baryons is simplistic. Clearly this problem needs to 
be tackled by high resolution hydrodynamical plus N-body 
simulations. We plan to undertake a detailed study of the 
impact of various astrophysics processes on this interaction 
in the future. In particular, this will affect our predictions 
for correlations involving the dynamics of galaxies (Tully- 
Fisher, Faber-Jackson and fundamental plane). 

• The modelling of dynamical friction relies on Chan- 
drasekhar's formula, which is correct on average but does 
not reflect the diversity of cases that appear in numerical 
simulations (Spring el et al. 2001| > . 

• Whilst we estimate the radial positions of galaxies in- 
side halos, the angular coordinates are not determined in 
this model. 

In particular, these points call attention to the need for more 
detailed/accurate modelling of what actually occurs to both 
gas and galaxies within cluster environments. 

10.1.2 Addressing the optical/ IR luminosity budget of 
galaxies 

We have introduced an efficient technique for modelling 
the distribution, composition and effects of dust in or- 
der to provide a full, panchromatic view of galaxy for- 
mation and evolution processes. For this purpose, we use 
the ingredients of the STARDUST model of chemical and 
spectrophotometric evolution of stellar populations (DGS). 
This model was specifically designed to be implemented 
within a semi-analytic model of galaxy formation and evo- 
lution ( Devriendt & Guidcrdoni 2000 1 . The star formation 
rate and cold gas metallicity histories are followed along the 
merging history trees of galaxies, and, once an IMF is cho- 
sen, the model gives theoretical spectra of the stellar pop- 
ulations. Given the galaxy radius and column density, the 
extinction and the IR/submm spectral energy distribution 



are computed without additional free parameters. Comput- 
ing the IR properties of galaxies is particularly important in 
the picture of hierarchical galaxy formation since observa- 
tions show that mergers radiate mostly in the IR. 

10.1.3 Studying resolution effects 

We provide some simple tests of the sensitivity of our re- 
sults to time and mass resolution. This enables us to put 
confident bounds on the regime where the results should be 
trusted. This is a very important step in the development of 
hybrid modelling, and resolution effects should be studied in 
more detail. The actual limit due to the baryonic mass of the 
lightest halos makes us incomplete fainter than Mb = —18.8 
at z — 0. This limitation is well understood and we simply 
have to appreciate that our predictive power is restricted 
to volume limited samples above this magnitude. However, 
the inability to resolve merging histories of small galaxies 
within this limit is a subtler problem, and requires more 
work to deal with. For instance, we attribute the underes- 
timate of the metallicity of the intra cluster medium to the 
absence of dwarf galaxies. Morover, the results of section l9~^l 
suggest that we are only completely accurate down to galax- 
ies around three-quarters of a magnitude brighter than this 
M B = -18.8 limit. 

10.1.4 Mimicking selection effects 

Our models produce a wide diversity of galaxies, but obser- 
vational data only probe the emerging part of this iceberg. 
We have stressed several times in this paper the need to 
understand observational selection effects when comparing 
models to data. Examples include the Malmquist bias in us- 
ing the APM catalogue to measure the relative abundance of 
different galaxy types, and the surface brightness selection 
that is generally present in all astrophysical datasets. We 
have roughly tested the influence of these selection effects 
on the comparison of the theoretical and observational sam- 
ples in these two cases. As a result, we manage to imitate 
the datasets in a plausible way. This deserves more work, 
and will be addressed in forthcoming papers. 

10.2 Main results 

While we have the flexibility to alter the free parameters 
governing the star formation rate, feedback efficiency, etc, 
we find that our observationally motivated initial choice of 
reference parameters is well able to match a number of ob- 
servational constraints: 

• We find a wide variety of galaxy properties, among 
which the Milky Way appears to have 'normal' properties for 
its circular velocity, in terms of baryonic mass, ISM metal- 
licity, star formation rate and size. 

• The luminosity functions in the optical and IR are well 
reproduced by the model, above the resolution limit. 

• The optical/IR luminosity budget in z — galaxies 
is well reproduced by the model, as it appears from the 
fit of the luminosity functions around at various opti- 
cal/IR wavelengths, and from the existence of an increasing 
IR/optical luminosity ratio with IR luminosity as observed 
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in the IRAS data. This is particularly interesting as the opti- 
cal/IR properties of galaxies are closely linked to the process 
of hierarchical galaxy formation, since mergers radiate most 
of their luminosity in the IR. These first results are thus 
very encouraging. 

• Because the definition of morphological types in usual 
catalogues of bright galaxies lies on a mix of 'objective' and 
somewhat 'subjective' classification criteria, it is not easy 
to compare the predictions of this kind of model with data. 
With our simple choice of the bulge-to-disc luminosity ratio 
as a tracer of the type, we fairly reproduce the type fractions 
in a volume-limited sample. We also reproduce the overall 
behaviour of the colour-type correlation, and find that the 
B — V colour monotonically reddens from the late to the 
early types. However we predict too slow an evolution when 
going from a morphological type to another, and therefore 
B — V colours which are too blue for the SOs/ellipticals. We 
are also unable to model the transition of the low-luminosity 
spirals to the 'dwarfs', that is expected to occur around 
Mb = —17, which is well below our resolution limit. 

• The sizes of high surface-brightness discs are in the 
same range as the observational data, but we predict slightly 
too many large bright discs than are observed. We recall 
that the sizes are computed under the assumption of spe- 
cific angular momentum conservation during the dissipative 
collapse of baryons, and in the absence of dynamical effects 
of the feedback mechanism on the structure of the galax- 
ies. Given the crudeness of the recipes which are used, we 
consider that this fit is still quite satisfying. 

• We obtain an arguably decent match of fundamental 
structural properties: the Tully-Fisher relation for the spi- 
ral galaxies, and the Faber-Jackson relation for the ellipti- 
cals. We compute the dispersion in these relations, and find 
good agreement for the TF of cD galaxies, but the scatter 
is too high and the slope of the relation itself too shallow 
for cluster satellites and field galaxies. The FJ scatter and 
slope for cluster early-type galaxies are in better agreement 
with those obtained from the data. 

• Finally, we obtain a fair fit to the fundamental plane in 
two wavebands, reproducing the observed parameters and 
the scatter. This is the first time semi-analytic models have 
been applied to look at this important observational test. 

Thus the impression that emerges from this first set of 
results confirms the overall robustness of the main statistical 
features predicted for local galaxies, with respect to previ- 
ous works involving semi-analytic models or the hybrid ap- 
proach, although the details of the current implementation 
in GALICS are fully new and original. Hierarchical galaxy for- 
mation does provide us with a natural explanation of many 
patterns of our local universe, at least qualitatively : the 
rough percentage of each morphological type along the Hub- 
ble sequence, the fact that early-type galaxies are redder on 
an average, the correlation between luminosity and a char- 
acteristic internal velocity, the existence of an upper limit 
for the luminosities of galaxies, and the overall luminosity 
budget between the optical and the IR. These patterns ap- 
pear very naturally from the interplay of a limited number 
of physical processes. 

The quantitative fit to local galaxies is generally satis- 
factory, given the crudeness of the assumptions, but some 
failures remain. More interesting perhaps is the inner con- 



sistency of these successes and failures. The massive galaxies 
seem to rotate too fast and be too blue, with a lot of ellipti- 
cals still having discs. This can be attributed to an improper 
modelling of cooling and/or to possible changes in the IMF, 
but it can be anticipated that a modification of the recipes 
will have simultaneous consequences on these three features. 
More work is needed from our reference model to clarify this 
point. 

The well-known influence of selection effects on the ob- 
served properties of galaxies receives a new emphasis from 
hierarchical galaxy formation. This theory generically pre- 
dicts that galaxies scatter within a wide range of proper- 
ties, and that observations always unveil only the tip of the 
iceberg. Such a pattern considerably complicates the com- 
parison of predictions with published data for which the 
exact observational selection criteria are not always fully 
documented. Large galaxy samples produced by extensive 
surveys, with a careful control of selection effects, should 
provide us with a better basis for the comparison with the 
predictions of the models. 

Other papers in the series will continue to detail the 
successes and failures of the current model at higher red- 
shifts, so that a consistent view on the evolution of most of 
the patterns can be obtained. 



11 CONCLUSIONS 

The semi-analytic approach has often been criticized for the 
number of free parameters it implies, and the differences 
between the various models. The impression after this fully 
original work here, when it is compared eg. to KCDW, is 
the overall robustness of many results, despite differences in 
the details of the recipes that are used to model the astro- 
physical processes, and to handle halos and galaxies. 

After a brief description of the iV-body simulation, and 
the way we build halo merging history trees, we have detailed 
the semi-analytic recipes that are used to model the baryons. 
We have then shown a series of predictions for galaxy proper- 
ties in the local Universe coming from a reference model. The 
model produces a wide range of properties of the local galaxy 
population, and gives reasonable fits to the optical/IR lumi- 
nosity functions, the correlation of the IR/optical luminosity 
ratio with IR luminosity, the morphological fractions, and 
the colour-morphology relation. Satisfactory fits of the disc 
sizes, the Tully-Fisher and Faber-Jackson relations, and the 
fundamental plane, are also obtained. We have tried to test 
the influence of various selection effects on the comparison 
of predictions to data. We have also shown that the influence 
of time and mass resolution on the results can be estimated, 
and should not affect our conclusions. 

One of the ambitions of the GALICS research programme 
is to make predictions for a full set of observational data, at 
various redshifts. Clearly such a programme will produce a 
host of results that will be covered by several papers in this 
series. Consequently, it is too early to form a full judgement 
of the ability of the model to capture the main statistical 
properties of galaxies, and their time evolution. In this first 
paper, we have presented our GALICS hybrid model for the 
study of galaxy formation. Forthcoming papers will comple- 
ment the picture given here. In paper II, we will study the 
influence of the cosmological and astrophysical parameters 
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on our results, and we will predict the redshift evolution of a 
set of statistical galaxy properties. Paper III will be devoted 
to Lyman-break galaxies. Paper IV will present multiwave- 
length faint counts and angular correlation functions, as well 
as the construction of mock images that are the most funda- 
mental basis for the study of selection effects. In paper V, we 
will study the 3D correlation function, and examine the be- 
haviour of the bias parameter for a variety of galaxy samples 
as a function of scale and time. After this global overview, 
subsequent papers will address more specific issues. 
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APPENDIX A: THE TREECODE 
Al The hierarchical approach 

When computing the gravitational potential acting on a 
particle, it is possible to ignore the details of the internal 
structure of a group of distant particles. Replacing several 
particle-particle interactions by one particle-distant group 
interaction is much more efficient than a particle-particle 
code, and one can still control the accuracy of the com- 
putations. This is the key simplification of the hierarchical 
approach. 

The first codes of this kind were written by |Appel| 
iPSTl GUiJ, |Jernigan (1985] , and |Porter (1985] , who 
were trying to follow the spatial distribution of mass as 
closely as possible. |Barnes &c Hut (1986} used a simpler 
algorithm, employing an hierarchical decomposition of 
space, allowing them to measure the accuracy of the force 
estimation. This code was vectorized by |Hernquist (1987| l, 
|Makino (1990} , and |Barnes (1990 1 and checked on vec- 
torial architectures like Cray X-MP. The adaptation 
for cosmological cases, especially the adding of peri- 
odic boundary conditions (using the Ewald summation 
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method) was performed by |Bouchet fc Hernquist (1 988), 
Hernquist, Bouchet, fc Suto (199lf | and 
Suginohara et al. (19910 . 

The tree-code uses a multi-polar expansion to compute 
the interaction of a particle and a group of particles. We can 
do this if the size of the group is small enough compared 
to the distance of the particle to this group: a variety of 
different criteria are in practice used for that condition. 

This concept is implemented by a hierarchical division 
of space, usually called a tree, which, in our case is a binary 
tree. The tree is constructed by starting from the entire box 
and dividing the space along one direction (say x first) in 
two equal cells. Each of those cells is divided along another 
direction (y), and this process (division along x, y, z, al- 
ternately) continues until each cell has one or zero particles 
inside. 

The computation of the force for one particle starts from 
the 'root' (that is the entire simulation box) of the tree. The 
selection criterion depending on the size of the cell as com- 
pared to the particle-cell distance is applied, to determine 
whether we can replace the cell by its centre of mass. If so, 
we add to the 'interaction list' the term corresponding to 
this interaction. If not, the same check is done, recursively, 
on sub-cells, until each term of the force is computed. 

This is the basic idea for the 'sequential' version of the 
code. We now describe the implementation of our parallel 
approach. 



A2 Parallelization 

Our code is written specifically for a T3D or T3E Cray ar- 
chitecture, which has a distributed memory type. Fast Cray 
libraries like ShMem allow us to use shared memory rou- 
tines, like remote reads or remote writes, which are rather 
faster than the usual message passing routines of PVM or 
even MPI. The code could be easily translated into MPI-2 
routines. 

At each timestep, a global sort of particles is performed 
on all processors. This global sort is based on four keys, 
which are defined by interleaving the binary digits of the 
x-, y-, and ^-positions of particles, in such a way that two 
physically close particles are, most probably, on the same 
processor, except for particles at the borders. This sorting 
operation thus amounts to a 'pre-building' of the tree, and 
also helps to decrease communications between processors 
during the force computation. The tree is then built, simul- 
taneously on all processors. The binary tree is a shared one, 
with root on processor 0. To decrease communications, cells 
that contain particles on the same processor will also be on 
this processor, most of the time; that is to say, going down 
along branches of the tree, we find that cells tend to be on 
the same processors as the particles they 'control'. 

The code also takes into account periodic boundary con- 
ditions by using the Ewald summation method, and copying 
the grid of periodic replicas onto each processor. 

The force computation is finally done in a parallel way 
using the shared memory library of the Cray T3E. Once 
forces are computed, we can derive accelerations, and then 
update particle velocities and positions, using a standard 
integration scheme. 

We describe a few finer details of the code below. 



• Sharing of particles on nodes. First of all, to bal- 
ance load, each processor is assigned approximately the same 
number of particles. More precisely, if -/Vp ar t is the total num- 
ber of particles, and JVprocs the number of processors, we 
perform the Euclidean division: 

Apart = qN plocs + r (Al) 

Each processor has at least q particles, and the first r pro- 
cessors have an additional one. 

• Sorting particles. Each particle is then assigned a 
key, composed of four integer numbers. Those keys are based 
on the idea of the bijection between R 3 and R. For com- 
puters, this is actually a bijection between N 3 and iV, be- 
cause of the limited precision of real numbers. The build- 
ing and sorting of particles is essentially based on Can- 
tor's construction of this bijection, which we briefly review 
here. We obtain the first three integers by interleaving the 
binary digits of the three spatial coordinates of the par- 
ticles, where all the coordinates are normalized between 
and 1. Let those coordinates be written as, for exam- 
ple: X = O.XlX 2 X3X4.X 5 X6X 7 X 8 X 9 , y = 0.2/1J/2J/3J/42/52/62/72/82/9, 

z = 0.ziZ2Z3Z4Z5ZeZ7Z8Zg. The first three integer keys will 
be: ki = x 1 y 1 z 1 x 2 y2Z2X3y3Z3, k 2 = X4y4Z4X 5 y 5 z 5 x 6 y6Z6, 
k3 = X7y7Z7Xsy&zsXgygZg. The fourth number used is just 
the original particle number, such that in the rare event that 
two particles have the same position, they will not have the 
same keys. 

The keys are used to sort particles on the processors in 
increasing lexicographic order. This is a global sort: after 
sorting, particles on processor number n + 1 will all have 
keys greater than those of the particles on processor n. 

• Building the tree. A global binary tree shared on 
all processors is then built. The entire box is the root of the 
tree. We then divide recursively all cells of the tree alter- 
nately along the (x, y, z) directions, until each cell is either 
empty or contains only one particle. As the tree is shared on 
all nodes, the top of the tree is divided on all processors. So, 
the beginning of the tree-building process requires a certain 
amount of communication between them. But, as the tree 
building proceeds, most of the sub-cells of a cell will be on 
the same processor (as the one controlling the cell), and this 
will finally decrease considerably inter-processor communi- 
cation. 

• Tree properties. We compute the mass and centre 
of mass positions of each cell. This is done by starting from 
the lowest cells in the tree, just above the particle level. A 
synchronisation of the computation is made between each 
level of the tree, because the results of the previous level are 
required to compute the mass and centre of mass position 
of this level. 

• Force computation. The force computation is done 
by starting from the top of the tree, going down recursively 
until each term of the force has been computed: each cell 
is either 'opened' or accepted in the interaction list, if the 
opening criterion is satisfied. 

• Particles update. A conventional leapfrog algorithm, 
taking into account the expansion of the Universe, is used 
to update positions and velocities of particles. 

• Optimisation. Several tricks have been used to opti- 
mize the code. Firstly we 'derecursified' most of the routines 
of the code. A copy of the top of the tree (the first N ~ 12 
levels), is also saved on each node, thus decreasing commu- 
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nications in the tree descent. To decrease cache misses, we 
also modify the tree, so that the right pointer previously 
used for the 'right child' is instead used, in the force com- 
putation, as a 'sibling pointer': this modification goes with 
the derecursification of the tree descent. 



APPENDIX B: HALO POTENTIAL ENERGIES 

As mentioned in section 12.31 the potential energy of 
large dark matter halos is computationally expensive us- 
ing a simple sum over pairs. We compare the measured 
potential with that from applying the formula used by 
|Steinmetz fc Bartelmann (1995| > for the potential energy of 
an ellipsoid, 

W = -\GM 2 R ¥ {a,b 2 ,c), (Bl) 

where _Rf is Carlson's form for the elliptic integral of the 
first kind (see eg. IPress et al. 19921 . 

R F (x,y,z) = - [ dt (B2) 

2 Jo y/(t + x)(t + y)(t + z) 

This procedure does not correct for the non-inverse square 
nature of the softened gravitational potential in the simula- 
tions, but, as shown in figure lrlTl there is a good correlation 
between the true potential energy and the energy estimated 
in this way. We note from this figure that the scatter in the 
relation is fairly small, of order five per cent. It will also 
be seen that the median depends on the number of parti- 
cles in the halo. It ranges from being a slight overestimate 
of the true potential energy at low halo mass, to a direct 
correspondence at around 25 particles, before flattening off 
towards large particle number with an underestimate of the 
magnitude of the potential energy of around fifteen per cent. 

In order to use an unbiased recipe, whilst still speed- 
ing up the computation, we apply this formula, to halos of 
more than one thousand particles, using the full double-sum 
method for smaller halos. The potential energies of these 
halos are then accurate only to five per cent. Since the spin 
parameter depends on the square root of the total energy, 
this additional scatter will have little effect on the distribu- 
tion of spin parameters. 

It should be noted that this cut-off means that the en- 
ergies are computed properly (ie. without this approxima- 
tion) for the majority of halos (given the mass distribution 
of Fig. , and for almost all the halos that were found to 
have positive total energy (since their mass function has a 
steeper slope). 




Figure Bl. Relationship between the true, measured potential 
energy from pair counting and that found by applying equa- 
tion lBll to the mass distribution. The faint dots show the relation 
for each halo, while the solid line is the median and the dotted 
lines the scatter at plus and minus one standard deviation. 



